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The  wastewater  chlorination  process  in  fixed-size  open-channel  reactors  is  a 
challenging  control  problem  because  of  the  presence  of  a variable  dead  time  caused  by 
changes  in  the  inlet  flow  rate.  The  approach  proposed  here  involves  the  utilization  of  a 
moving  weir  as  an  actuator  in  a novel  residence  time  control  scheme,  designed  to  ensure 
that  fluid  passing  through  the  reactor  experiences  a constant  residence  time  independent 
of  flow  rate.  The  dynamic  weir  control  loop  is  conceived  to  supplement  existing  control 
structures  that  are  charged  with  manipulating  the  dosage  of  chlorine.  Simulations  studies 
show  that  both  feedforward  and  feedback  manifestations  of  the  residence  time  controller 
are  highly  effective  at  improving  chlorine  dosage  control  and  are  highly  robust  with 
respect  to  weir-flow  modeling  and  measurement  errors.  In  fact,  these  controllers  have  the 
potential  to  improve  the  control  performance  dramatically  while  saving  as  much  as  45% 
of  the  chlorine  used  for  treatment. 
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An  automatic  control  system  described  as  run-to-run  is  one  in  which  a discrete 
controller  is  applied  to  an  inherently  discrete  process,  and  where  the  performance 
measurements  are  made  only  after  all  processing  is  completed.  A complete  theoretical 
analysis  is  developed  for  the  case  of  a static  plant  under  proportional-integral  control, 
resulting  in  quantitative  descriptions  of  the  nominal  and  robust  stability  regions,  and 
recommended  tuning  heuristics.  The  analysis  for  the  case  of  a Smith  Predictor  design 
shows  that  very  small  modeling  errors  in  the  estimated  measurement  delay  may  cause 
closed-loop  instability.  Therefore,  the  use  of  a Smith  Predictor  for  run-to-run  systems 
suffering  from  uncertainty  in  the  delay  is  not  recommended. 

A weak-acid/strong-base  neutralization  control  system  utilizing  two  alternative 
pH  sensors  is  modeled  and  analyzed.  The  first  sensor  is  a traditional  pH  electrode  whose 
output  signal  is  proportional  to  pH,  and  the  second  sensor  is  a fiber-optic  optrode  that 
produces  a signal  that  is  proportional  to  hydrogen  ion  concentration.  Systematic 
simulation  studies  show  that  no  inherent  performance  improvement  is  observed  when 
using  the  optrode  instead  of  the  electrode  sensor  in  the  feedback  loop. 


IX 


CHAPTER  1 
INTRODUCTION 


There  are  three  types  of  process  control  studies  addressed  in  this  dissertation, 
namely,  a design  project,  a theoretical  project,  and  an  experimental  project.  The  design 
project  examines  wastewater  chlorination  control,  a distributed-parameter  system  that  is 
typically  difficult  to  control  with  standard  strategies  due  to  the  presence  of  high 
variations  in  flow  rate.  It  is  shown  that  a novel  design  modification  based  on  an 
automatically  adjusted  moving  weir  can  lead  to  dramatically  improved  performance.  In 
the  theoretical  study,  single-input  single-output  run-to-run  unit  operations  are 
investigated  for  the  case  of  a static  plant.  The  run-to-run  controller  form  is  contrasted 
with  a standard  discrete  controller  form,  exposing  subtle  differences  in  behavioral 
characteristics  under  closed-loop  control.  Finally,  the  experimental  project  seeks  to 
improve  the  performance  of  pH  control  systems  that  use  traditional  pH  electrodes. 

Known  to  be  highly  non-linear,  the  electrode  may  have  limited  performance  potential  as  a 
result.  A novel  linear  optical  sensor  is  introduced  as  an  alternative. 

Chapter  2 presents  a design  study  involving  a standard  wastewater  chlorination 
basin.  It  is  shown  that  traditional  control  schemes  have  difficulty  maintaining  the  outlet 
concentration  at  a specified  value.  The  difficulty  is  caused  by  the  effect  that  the  varying 
flow  rate  of  wastewater  has  on  the  deadtime  (residence  time)  of  the  system.  Because  the 
volume  of  the  reactor  is  constant,  residence  time  is  inversely  proportional  to  flow  rate, 
and  the  latter  can  vary  by  as  much  as  an  order  of  magnitude.  Specifically,  the  variable 
residence  time  causes  poor  closed-loop  performance.  The  system  is  modified  by 
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installation  of  a moving  weir,  which  allows  a controller  to  manipulate  the  volume  of  the 
reaction  chamber.  When  matched  to  the  flow  rate  of  wastewater,  the  residence  time  can 
be  held  approximately  constant. 

To  validate  the  effectiveness  of  the  design  modification,  a comprehensive 
simulation  of  the  system  is  used.  The  first  step  of  the  process  is  the  development  of  a 
realistic  representation  of  the  flow  through  a basin.  Next,  traditional  control  schemes  are 
considered,  and  their  performance  is  evaluated.  Finally,  the  basin  model  basin  is 
modified  to  allow  for  the  automatic  adjustment  of  the  volume  through  use  of  a moving 
weir,  and  controller  performance  is  again  tested.  It  is  at  this  point  that  the  merit  of  the 
design  modification  is  clear,  as  it  effectively  rejects  the  given  disturbance.  Because  the 
total  chlorine  consumption  is  highly  influenced  by  residence  time,  the  movable  weir 
facilitates  significant  reduction  in  chlorine  usage  while  attaining  highly  improved  outlet 
concentration  control  by  the  existing  control  scheme. 

A theoretical  project  that  is  the  subject  of  Chapter  3 addresses  the  analysis  of  run- 
to-run  systems  where  significant  measurement  delay  is  observed.  A run-to-run  system  is 
one  in  which  a unit  operation  processes  batches  of  material  in  succession,  with  the 
constraint  that  each  batch  must  be  completed  before  the  relevant  end-product  properties 
may  be  determined.  Examples  include  a wide  range  of  processes,  such  as  industrial 
peanut  roasting  and  compact  disc  sputtering  operations.  In  each  case,  the  unit  operation 
that  imparts  some  physical  characteristic  on  each  batch  is  only  one  of  many  in 
succession.  Consequently,  there  is  significant  delay  observed  by  the  time  said 
characteristic  can  be  measured,  and  only  then  can  adjustments  be  made  to  the  system.  An 
important  finding  of  the  work  is  that,  although  run-to-run  systems  have  characteristics  in 
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common  with  traditional  discrete  systems,  they  behave  markedly  different  under  closed- 
loop  conditions,  and  care  must  be  taken  by  the  control  designer  to  guarantee  that  the 
system  will  remain  stable  under  realistic  conditions. 

The  mathematical  description  is  developed  using  standard  techniques,  primarily, 
the  Jury  Array,  which  provides  necessary  and  sufficient  conditions  to  guarantee  that  the 
roots  of  a polynomial  are  inside  the  unit  circle.  Discrete  systems  are  necessarily  unstable 
if  at  least  one  (uncancelled)  pole  lies  outside  the  unit  circle,  making  the  Jury  Array  an 
appropriate  stability  test  for  this  type  of  system.  Next,  because  significant  measurement 
delay  is  typical  of  many  run-to-run  systems,  the  performance  of  a deadtime  compensator 
is  analyzed.  It  is  shown  that  such  a compensator  is  ineffective,  primarily  because  a run- 
to-run  system  is  fundamentally  different  from  the  traditional  discrete  system  for  which 
compensators  are  designed. 

Chapter  4 presents  a thorough  investigation  of  the  properties  of  linear  optical  pH 
sensors.  In  contrast,  traditional  pH  electrodes  are  logarithmic  with  respect  to  the 
measured  variable,  namely,  the  hydrogen  ion  concentration.  Consequently,  the  two  types 
of  sensors  behave  differently  in  closed-loop  feedback  control  loops.  Generally,  the  linear 
sensor  is  expected  to  perform  better  than  a non-linear  one,  and  the  primary  result  of  the 
study  is  to  quantify  the  design  specifications  in  terms  of  how  much  the  physical 
properties  of  the  optical  sensor  can  degrade  while  maintaining  equivalent  or  better 
performance  than  the  electrode.  Similarly,  it  is  also  shown  how  much  improvement  can 
be  expected  by  using  the  optical  sensor  assuming  a realistic  set  of  physical 


characteristics. 
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Both  computer  simulations  and  bench-scale  experiments  are  used  in  the  pH 
sensing  study.  Initially,  the  important  chemical  reactions  that  take  place  in  an  acid-base 
neutralization  process  are  modeled  and  simulated.  Next,  the  bench-scale  pH  control 
system  is  utilized  for  the  purpose  of  verifying  the  validity  of  the  simulation  behavior. 
Simulations  are  used  for  the  final  phase  of  the  project,  which  quantifies  the  perfonnance 
of  the  optical  sensor  with  respect  to  standard  electrodes. 


CHAPTER  2 

WASTEWATER  CHLORINATION  USING  A NOVEL  CONTROL  SCHEME 

2.1  Introduction 

A major  requirement  of  wastewater  treatment  plants  is  that  the  treated  water  must 
be  effectively  disinfected  before  it  is  discharged.  Although  this  can  be  done  in  several 
ways,  the  most  common  approach  is  through  chemical  reaction  with  chlorine  (Metcalf  & 
Eddy,  1979).  Chlorine  is  consumed  in  the  process,  but  regulations  stipulate  that  a 
residual  chlorine  concentration  must  be  present  at  or  above  a minimum  level  upon 
discharge  from  the  reactor.  The  restriction  ensures  that  enough  chlorine  is  present  in  the 
water  for  a long  enough  time  to  guarantee  satisfactory  cleansing  (Hammer  and  Hammer, 
1996).  However,  excessively  high  chlorine  levels  are  problematic  for  three  reasons. 

First,  the  purchase  of  chlorine  is  a significant  operating  cost.  Second,  unnecessary 
addition  of  chlorine  means  that  greater  volumes  must  be  stored  on-site,  which  can  be  both 
costly  and  dangerous.  Third,  under  high  chlorine  concentrations  the  formation  of 
trihalomethane  (THM)  contaminants  is  favored  (Jolley  et  al.,  1990),  and  accordingly,  the 
process  must  satisfy  regulations  restricting  their  maximum  level.  It  is  therefore  desirable 
to  maintain  the  outlet  chlorine  concentration  as  closely  as  possible  to  a predetermined 
value. 

The  goal  of  this  research  is  to  develop  an  effective,  automatic,  chlorine-dosage 
control  scheme  for  a wastewater  treatment  plant.  With  relatively  little  storage  capacity, 
many  plants  must  process  incoming  water  as  it  arrives,  and  are  therefore  subjected  to  the 
daily  cyclical  nature  of  water  usage  (Hammer  and  Hammer,  1996).  Because  a reactor’s 
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size  is  fixed,  the  varying  rate  of  the  flowing  water  poses  the  most  significant  technical 
challenge  to  automatic  control — it  introduces  a time-varying  system  delay,  which 
significantly  degrades  treated  wastewater  quality  (Smith,  1971).  In  addition,  other 
disturbances  such  as  intense  solar  radiation  (causing  increases  in  volatilization  of 
chlorine,  evaporation  of  water,  and  chemical  reaction  rates)  and  rainfall  (causing  dilution) 
can  affect  the  dynamics  of  these  outdoor  reactors. 

A number  of  feedback  control  designs  can  be  considered  as  candidates  to  ensure 
optimal  outlet  chlorine  concentration  levels.  The  least  complex  strategy  is  simple 
proportional-integral  (PI)  feedback,  but  this  approach  is  usually  inadequate  for  processes 
with  large  input  disturbances  and  deadtime  (Seborg,  Edgar,  and  Mellichamp,  1989).  At 
the  price  of  increased  complexity,  a cascade  regime  helps  minimize  the  effects  of  input 
disturbances.  Furthermore,  adding  a Smith  Predictor  (which  requires  an  accurate  model 
of  the  reactor)  can  counter  the  adverse  effects  of  deadtime,  provided  the  deadtime  has 
limited  variability  (Ogunnaike  and  Ray,  1994).  While  many  of  these  attributes  are 
promising,  it  is  unlikely  that  they  are  fully  adequate,  even  in  combination,  to  control  a 
system  with  such  a highly  varying  deadtime.  It  is  proposed  that  a process  modification 
may  prove  effective  in  aiding  the  control  of  residual  chlorine  concentration.  If  the  end  of 
the  reactor  is  fitted  with  a dynamic  weir,  the  volume  of  the  reactor  could  be  varied;  this 
would  allow  the  operator  to  manipulate  the  residence  time  of  the  fluid  in  the  reactor  and, 
hence,  potentially  ensure  that  the  time-delay  of  the  system  remains  approximately 
constant. 

An  accurate  model  of  a wastewater  treatment  plant  can  be  used  to  quantify  the 
relative  effectiveness  of  each  candidate  control  approach.  The  sections  contained  in  this 
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chapter  identify  issues  pertinent  to  the  formulation  of  a viable  model,  as  well  as  the 
results  of  applying  various  control  schemes.  More  specifically,  Section  2.2  describes  all 
the  major  aspects  of  reactor  modeling.  Critical  parameters  of  the  derived  reactor  model 
are  then  identified  in  Section  2.3  using  experimental  plant  data.  Next,  the  design  and 
performance  assessment  is  presented  in  Section  2.4  for  several  possible  control  schemes. 
Then,  Section  2.5  details  the  potential  of  a novel  process  modification — a dynamic 
weir — as  an  alternative  design,  and  finally,  conclusions  for  the  most  effective  strategy  to 
control  residual  chlorine  concentration  in  wastewater  are  outlined  in  Section  2.6. 

2.2  Process  Model  of  a Typical  Wastewater  Disinfection  Plant 
An  abstract  representation  of  a disinfection  system  is  depicted  in  Figure  2-1 . The 
rectangle  at  the  left  of  the  figure  represents  the  dosing  area,  where  chlorine  is  mixed  into 
the  system.  The  chlorine  reacts  with  ammonia  nitrogen  (and  other  trace  species)  present 
in  the  incoming  water,  consuming  a portion  of  the  chlorine  very  quickly.  The  amount 
consumed  by  this  process  is  referred  to  as  the  ammonia  nitrogen  chlorine  demand  (the 
effect  of  the  trace  species  is  included  in  this  term  in  this  study  for  brevity).  The  dosing 
area  is  designed  to  eliminate  ammonia  nitrogen  from  the  water  and  to  promote  an  outlet 
flow  with  a uniform  free-plus-combined  (FPC)  residual  chlorine  concentration  (Metcalf 
& Eddy,  1979).  The  outlet  flow  immediately  enters  the  next  stage,  or  contact  area,  where 
the  bulk  of  disinfection  takes  place.  There  are  four  important  values  associated  with 

determining  the  chlorine  concentration  in  the  dosage  area:  namely,  Q(t),  the  influent 

3 

volumetric  flowrate  measured  in  cubic  meters  per  second  (m  Is);  5(t),  the  ammonia 
nitrogen  chlorine  demand  of  the  incoming  water  measured  in  parts  per  million  (ppm) 

3 

(where  ppm  is  g/m  );  V(]a,  the  volume  of  liquid  in  the  dosage  area  measured  in  cubic 
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meters  (m  );  and  u(t),  the  chlorine  dosage  input  rate  measured  in  grams  per  second  (g/s). 

Water  leaving  the  dosage  area  has  a concentration  Co(t)  measured  in  parts  per  million 
(ppm).  The  dosing  area  is  distinguished  from  the  contact  area  at  a point  where  uniform 
chlorine  concentration  Co(t)  and  plug-flow  can  be  guaranteed.  Physically,  the  dosing  and 

contact  areas  are  parts  of  the  same  basin;  there  may  or  may  not  be  some  partition  (such  as 
a gate)  between  the  two  areas,  depending  upon  the  design  of  the  plant. 

The  rectangle  on  the  right  of  the  figure  represents  the  contact  area,  where  the  FPC 
residual  chlorine  disinfects  the  water.  Disinfecting  chlorine  species  are  slowly  consumed 
by  chemical  reactions  during  the  time  the  liquid  resides  in  the  contact  area,  as  well  as  lost 
through  volatilization  to  the  surrounding  air.  In  contrast  to  the  dosage  area,  minimal 
mixing  is  desired  in  the  contact  area,  thus  ensuring  that  all  fluid  that  passes  through  it  is 
retained  for  a length  of  time  sufficient  for  disinfection.  The  contact  area  features  two 

relevant  variables:  namely,  Q(t),  as  defined  before;  and  Vca,  the  volume  of  the  liquid  in 

3 

the  contact  area  measured  in  cubic  meters  (m  ).  In  addition,  a z-coordinate  is  introduced 

to  represent  the  longitudinal  direction  of  fluid  travel.  This  facilitates  reporting  the  value 
of  variables  at  various  spacial  locations.  Water  passing  through  the  contact  area  enters 

with  a concentration  Co(t),  moves  in  the  direction  of  increasing  z while  chemical 
reactions  and  volatilization  take  place,  and  exits  the  contact  basin  at  the  coordinate  point 
z = L with  a concentration  Cl(0.  The  reactor  is  normally  operated  at  constant  volume  by 

matching  the  inlet  and  outlet  flowrates.  Hence,  the  liquid  volumes  and  Vca  are 
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constants.  The  total  volume  of  the  liquid  in  the  contact  area  is  given  by  Vca  = AL,  where 
A is  the  cross-sectional  area  of  the  liquid  volume. 


Co  cL 


Q(t) 

u(t) 

w 

Dosage  Area 

Contact  Area 

m 

Vda 

Vca 

► 

Figure  2-1.  Schematic  of  a wastewater  treatment  plant, 
operated  at  constant  volume. 

Water  disinfection  plants  implement  the  abstract  design  of  Figure  2-1  in  various 
forms.  For  example,  space  considerations  often  dictate  that  the  contact  basin  is 
serpentine  with  a baffle-style  arrangement.  However,  as  long  as  uniform,  plug-flow  is 
preserved,  such  variations  are  conceptually  equivalent  to  the  abstract  linear  model 
proposed.  The  Kanapaha  water  treatment  facility  in  Gainesville,  Florida,  is  the  basis  for 
the  process  model  adopted  in  this  study,  and  Table  2-1  lists  various  physical  quantities  of 
the  plant,  used  for  both  modeling  and  model  validation  purposes. 


Table  2-1.  Physical  quantities  that 
describe  the  reactor  at  Kanapaha  water 
treatment  facility  in  Gainesville,  Florida. 


V 

2880  m3 

Vda 

57.6  m3  (0.02  V)  f 

Vca 

2822.4  m3  (0.98  V) 

j L 

1 20  m (average) 

A 

2 

24  m (average) 
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In  addition,  Table  2-2  lists  various  typical  operating  conditions  for  the  same 
facility.  The  first  three  quantities  are  the  typical  chlorine  concentration  at  the  entrance  of 
the  reactor  u',  the  dosing-area/contact-area  boundary  concentration  C0,  and  the  exit 
concentration  of  the  reactor  CL , measured  in  parts  per  million  (ppm).  The  final  quantity 
is  the  typical  flowrate  Q through  the  plant  measured  in  cubic  meters  per  second  (m  /s). 

Note  that  the  dosing  rate  u (measured  in  g/s)  is  fundamentally  different  from  the  dosage 
concentration  u'  (measured  in  ppm);  this  is  an  important  distinction  that  is  discussed  in 
Section  2.2.1. 


Table  2-2.  Typical  operating  conditions 
for  the  reactor  at  Kanapaha  water 
treatment  facility  in  Gainesville,  Florida. 


u' 

8.5  ppm 

Co 

5.5  ppm 

CL 

2.5  ppm 

j 

3 

Q 

0.4  m /s 

5 

3.0  ppm 

2.2.1  Kinetic  Model  of  Chlorine  Reaction  in  a Wastewater  Treatment  Plant 

Using  chlorine  to  eliminate  ammonia  and  disinfect  water  is  complicated  in  terms 
of  chemistry.  There  is  an  initial  series  of  reactions  involving  chlorine,  ammonia  nitrogen, 
and  their  products  in  the  dosing  area.  Some  of  those  products  then  act  as  disinfecting 
agents,  and  are  slowly  consumed  as  the  water  moves  down  the  length  of  the  contact 
basin.  The  ammonia  reactions  are  described  by  the  following  Morris-Wei  model  (Morris 
and  Wei,  1969): 


NH3  + HOC1  ->  nh2ci  + h2o 


(2-1) 


NH2C1  + HOC1  ->  nhci2  + h2o 

(2-2) 

nhci2  + HOC1  ->  nci3  + h2o 

(2-3) 

nci3  + H20  ->  NHC12  + HOC1 

(2-4) 

NHC12  +H20->N0H+2H++2Cr 

(2-5) 

noh+nh2ci  -»  n2  +h2o  + h+  +cr 

(2-6) 

NOH  + NHCl2  ->N2+HOCl  + H++Cr 

(2-7) 

NOH  + 2HOC1  ->  N03  + 3H+  + 2CF 

(2-8) 

These  reactions  proceed  rapidly  to  steady-state,  with  some  reactions  going  to 
completion,  and  others  reaching  equilibrium.  It  is  unclear  why  some  reverse  reactions 
(such  as  (2-4))  are  not  noted  in  the  model,  but  it  is  certain  that  at  least  nine  different 
chemical  species  are  involved.  In  addition,  the  reaction-rate  constants  are  also  either 
temperature  or  pH  dependent,  or  both.  For  example,  pH  and  temperature  strongly 
influence  the  point  of  equilibrium  for  the  chloramine  reactions  (2-1)  to  (2-4)  (Hammer 
and  Hammer,  1996).  While  the  Morris-Wei  model  is  a well-founded  proposal  of  what 
chemical  reactions  take  place  during  ammonia  destruction,  it  introduces  numerous 
differential  and  algebraic  relationships,  hence  complicating  the  theoretical  analysis  and 
increasing  the  computational  burden  significantly.  Furthermore,  concentrations  of  some 
of  the  nine  species  described  are  additionally  influenced  by  organic  nitrogen,  metal  ion, 
and  other  gas-forming  reactions  not  represented  in  the  Morris-Wei  model. 
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It  is  therefore  desirable  to  find  a simplified  though  realistic  model.  The 
knowledge  that  the  ammonia  is  almost  completely  consumed  very  rapidly  may  be 
exploited  in  deriving  a simplified  model.  One  possible  approach  is  to  treat  the  dosing 
area  as  a CSTR  within  which  the  ammonia  reactions  take  place.  The  resulting  model  is 
described  in  Section  2.2.2.  Although  the  dosing  area  is  not  always  physically  similar  to  a 
CSTR  (for  example,  it  can  include  a diffuser,  or  a line  source  of  chlorine  being  injected 
into  the  influent),  it  is  argued  that  the  relevant  behavior  of  the  dosing  area  is  adequately 
reflected  by  this  model,  with  the  added  benefit  that  the  model  is  conceptually  quite 
simple. 

The  disinfection  reactions  in  the  contact  area,  however,  are  not  intended  to 
proceed  to  completion  and  do  not  necessarily  reach  equilibrium.  The  regulations  that 
require  minimum  residual  chlorine  refer  specifically  to  free-plus-combined  (FPC) 
chlorine,  which  are  the  four  species  that  promote  disinfection.  They  include 

monochloramine  (NH2CI),  dichloramine  (NHCI2),  hypochlorous  acid  (HOC1),  and 
hypochlorite  ion  (OC1  ).  The  four  are  collectively  referred  to  as  “chlorine”  with  regard  to 

the  contact  area,  and  are  generically  represented  by  the  symbol  Cl  . Because  disinfection 

chemistry  is  even  more  complicated  than  ammonia  destruction  chemistry,  it  is  again 
desirable  to  find  a simplified  description  of  chlorine  concentration.  It  is  proposed  that  the 
chlorine  consumption  in  the  contact  area  can  be  modeled  through  the  lumped  chemical 
reaction 


Cl  — > Products  (2-9) 

with  a phenomenological  first-order  overall  rate  of  reaction  (Fogler,  1992) 
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rrxn  krxnC  (2-10) 

where  krxn  is  the  overall  reaction  rate  constant  of  chlorine  measured  in  reciprocal  seconds 

(1/s),  and  C is  the  concentration  of  chlorine  measured  in  parts  per  million  (ppm).  In 
addition,  FPC  chlorine  disappears  from  the  reactor  through  volatilization  and  undesired 
chemical  reactions  dependent  upon  UV  radiation  exposure;  these  phenomena  occur  at 
rates  that  are  also  concentration-dependent.  Therefore,  it  is  plausible  to  account  for  these 
effects  within  the  lumped  parameter  approach  of  (2- 1 0).  It  is  initially  assumed  that  a 
first-order  model  is  appropriate,  since  most  chemical  species  that  are  present  in  any 
appreciable  quantity  over  the  operating  range  of  the  system  are  monatomic  in  chlorine 
(Snoeyink  and  Jenkins,  1980).  The  suitability  of  the  overall  kinetic  expression  (2-10)  is 
assessed  using  experimental  data  in  Section  3. 

2.2.2  Dosage  Area  Model 

A first-principles  model  for  the  dosage  area  is  based  upon  the  constant-volume 
continuously-stirred  tank  reactor  (CSTR)  formulation  (Ogunnaike  and  Ray,  1994) 

v'"  = °" (,)Cj  “ Q(t,c»(1) + rV*  (2-ii > 

dA^  = Q(t)ANm  (t)  - Q(t)AN(t)  + rAN Vd,  (2-12) 

3 

where  Q<j(t)  is  the  flowrate  of  the  manipulated  dosage  stream  (m  /s),  Cj  is  the  constant 

3 

chlorine  concentration  in  the  dosage  stream  (g/m  ),  r is  the  rate  of  reaction  of  chlorine 

3 

(g/s),  ANin  is  the  effective  incoming  ammonia-nitrogen  concentration  (g/m  ),  AN  is  the 

3 

effective  ammonia-nitrogen  concentration  exiting  the  reactor  (g/m  ),  and  rAN  is  the  rate 
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of  reaction  of  ammonia  nitrogen  (g/s).  The  term  “effective”  indicates  that  the  bulk  of  the 
reactant  is  ammonia  nitrogen,  with  the  remainder  comprised  of  trace  reactants  such  as 
metal  ions. 

A series  of  assumptions  and  observations  allows  the  two  equations  to  be 
combined  into  one.  First,  although  the  slow  reactions  begin  in  the  dosing  area,  the 
residence  time  is  low  enough  that  they  are  considered  negligible,  so  they  are  ignored. 
Second,  note  that  r = R(t)rAN  , where  R is  the  ratio  of  the  stoichiometric  coefficients  for 
the  effective  ammonia  and  chlorine  reactions  (after  appropriate  conversion  from  molar 
units  to  mass  units).  Note  that  this  term  is  time  dependent,  since  the  trace  reactant 
concentration,  pH,  and  temperature  may  cause  the  demand  for  chlorine  for  each  unit  of 
ammonia  to  fluctuate.  In  other  words,  every  gram  of  ammonia  reacts  quickly  to 
equilibrium  with  R grams  of  chlorine,  and  R is  dynamic.  Also,  AN  = 0,  because  it  has 
already  been  assumed  that  virtually  all  of  the  ammonia  is  consumed  in  the  CSTR.  Next, 
dAN  n . 

= 0 , since  it  is  assumed  that  the  ammonia-nitrogen  quickly  reaches  its  steady-state. 

Finally,  it  should  be  noted  that  the  assumption  r = 0 cannot  be  made,  because  although 
rc„  = kC0nAN‘"  and  AN  = 0,  the  value  of  k is  very  large,  and  therefore  prevents  the  term 

from  being  effectively  zero.  The  steady-state  and  zero  ammonia-nitrogen  concentration 
assumptions  reduce  (2-12)  to  rAN  Vda  = Q(t)AN(t) , and  then  using  r = RrAN  to  substitute 
into  (2-11)  yields 


Vda  = Qd  (t)Cd  - Q(t)C0  (t)  - Q(t)R(t)AN(t)  (2-13) 


Recognizing  that  both  R(t)  and  AN(t)  are  unmeasured  disturbances,  let 
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R(t)AN(t)  = 8(t)  (2-14) 

where  8(t)  is  referred  to  as  the  chlorine  demand , is  measured  in  ppm  (g/s),  and  defines 

how  many  ppm  of  chlorine  are  consumed  by  fast  reactions  per  unit  of  flow.  When  this 
demand  is  satisfied,  the  system  is  said  to  have  reached  breakpoint  (Bierck  et  al.,  1996). 
Therefore,  the  final  CSTR  model  is  cast  in  the  form 


Vd3  = Qd  (t)Cd  - Q(t)C0(t)  - Q(t)8(t)  (2-1 5) 

dt 

The  model  maintains  a high  level  of  simplicity  while  capturing  the  important 
characteristics  of  a dosing  area.  First,  it  reflects  that  two  streams  are  being  combined  to  a 
uniform  concentration.  In  addition,  it  captures  the  fact  that  some  chemical  reactions 
quickly  proceed  to  completion  and/or  equilibrium;  the  resulting  behavior  can  be  viewed 
as  a load  on  the  system.  Finally,  it  is  assumed  that  any  differences  in  the  dynamic 
behavior  of  a CSTR  and  an  actual  dosing  area  are  negligible  on  the  time  scale  of  the 
contact  area,  since  its  volume,  and  therefore  its  residence  time,  typically  comprises  only  a 
small  fraction  (perhaps  1 or  2 percent)  of  the  total  volume  of  the  plant. 

A standard  linearization  (Seborg,  Edgar,  and  Mellichamp,  1 989)  of  the  model  (2- 
15)  yields 


dC0(t) 

dt 


c0(t)+ 


Cd. 

Vda 


QdW- 


^(t)_(Co±S) 
Vda  Vda 


Q(t) 


where  the  deviation  variables 


(2-16) 


C0:=C0(t)-C0 

(2-17) 

Qd  ;=  Qd(t)~Qd 

(2-18) 

5 : = 5(t)-8 

(2-19) 

16 


Q:=Q(t)-Q  (2-20) 

are  defined  as  a function  of  the  steady  state  values  C0,  Qd,8,  and  Q that  make  the  right- 

hand  side  of  (2-16)  identically  zero.  The  Laplace  transform  of  (2-16)  yields  the  transfer- 
function  model 


co(s)  = T7-2—  Qd(s)- S(s)- 

Vnu  * V , 


4^-s  + l 


Q Q 

Note  that  since  the  manipulated  variable  is 


i^s  + 1 


C0  +5 


^s  + l 

Q 


Q(s) 


(2-21) 


then 


where 


and  (2-21)  becomes 


u(s):=Qd(s)Cd 


u(s):=  Qd(s)Cd 


u(s):=  u(s)-  u 


C0(s)=  v 

* i 


da 

Q 


u(s)- 


1 


s + 1 


^s  + l 

Q 


5(s)- 


C0+5 


Vda 

-^S  + 

Q 


Q(s) 


(2-22) 


(2-23) 


(2-24) 


(2-25) 


Values  of  the  constants  appearing  in  (2-25)  are  given  in  Table  2-2. 

It  is  common  practice  in  industry  to  report  the  chlorine  dosage  in  ppm  (as  in  Table 
2-2),  which  is  calculated  by  dividing  the  flowrate  Q(t)  into  the  dosing  rate  Qd(t)Cd . The 


goal  is  to  report  u(t),  Cp(t),  and  Cp(t)  all  in  units  of  ppm.  The  result  introduces  the 


alternate  dosage  variable 
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u'(t) 


•=  Qd(0Cd 
Q(t) 


(2-26) 


where  u'(t)  is  measured  in  ppm  (/.<?.,  g/m  ).  The  product  Qd(t)Cd  is  the  equivalent  of 

u(t),  the  mass  dosing  rate  of  chlorine,  and  (2-26)  remains  valid  provided  that 
Qd(0  <<c  Q(t),  which  is  the  case  at  all  typical  operating  conditions.  This  is  because, 
strictly  speaking,  the  chlorine  dosing  rate,  Qd(t)Cd , should  be  divided  by  the  sum  of  the 
flowrates,  Q(t)  + Qd  (t) , to  calculate  the  concentration;  it  is  assumed  that  the  flowrate 
Qd(t)  is  negligible  compared  to  Q(t) . However,  using  (2-26)  would  introduce  a bilinear 
term  to  (2-15)  and  (2-25),  which  would  result  in  a more  complex  model,  and  so  u'(t)  is 
used  throughout  the  study  strictly  for  reporting  purposes,  while  u(t)  remains  as  the  actual 
control  variable.  Of  course,  using  the  sum  of  the  flowrates  in  the  denominator  of  (2-26) 
would  be  even  more  non-linear,  and  is  also  therefore  eliminated  from  consideration  as  a 
choice  of  control  variable. 

2.2.3  Contact  Area  Model 

The  most  important  component  of  the  system  is  the  contact  basin,  and  its  model  is 
therefore  critical.  As  a distributed-parameter  system  (Seborg,  Edgar,  and  Mellichamp, 
1989),  the  contact  area  can  be  modeled  as  a plug  flow  reactor  (PFR)  with  chemical 
reaction  through  the  standard  equation 


ac  ac  ^ a2c 

¥ = -v(,)fe+r“"  + D^7 


(2-27) 


where  C - C(t,z)  is  the  chlorine  concentration  measured  in  parts  per  million  (ppm), 
v(t)  = Q(t)  / A is  the  linear  flow  velocity  measured  in  meters  per  second  (m/s),  t is  the 

time  measured  in  seconds  (s),  rmx  is  the  rate  of  chemical  reaction  measured  in  parts  per 
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million  per  second  (ppm/s),  z is  the  horizontal  length  coordinate  of  the  reactor  measured 
in  meters  (m),  and  D is  the  effective  dispersion  coefficient  measured  in  meters  squared 
per  second  (m2/s).  A complete  derivation  of  Equation  (2-27)  is  given  in  Appendix  A. 

The  boundary  condition  is  C(t,0)  = Co(t),  and  the  initial  condition  is  C(0,z)  = C°(t,z), 
where  C (t,z)  is  a known  initial  spatial  concentration  profile. 

Three  primary  modeling  issues  arise.  First,  the  rate  of  chemical  reaction  rmx  must 

be  adequately  described;  this  matter  is  investigated  in  Section  2.2.3.  Second,  Equation 
(2-27)  is  a partial  differential  equation  (PDE),  which  is  not  directly  useful  for  the 
purposes  of  control  design;  Sections  2.2.4  and  2.2.5  outline  one  approach  to  a solution. 
Third,  since  an  ideal  PFR  passes  any  input  proportionally  after  a delay,  and  the  delay  is 
defined  as  the  residence  time  0ca(t)  = Vca  / Q(t),  the  model  predicts  that  this  system 
features  a delay  that  is  inversely  proportional  to  the  flowrate.  Because  the  flowrate  is 
time-varying,  it  is  critical  to  select  a control  design  that  can  perform  effectively  even  if 
the  delay  changes.  The  performance  for  different  schemes  is  explored  in  Section  2.4. 

2.2.4  Simplified  Model  of  a Plug  Flow  Reactor 

Substituting  the  kinetic  rate  model  (2-10)  and  the  linear  velocity  expression  v(t)= 
Q(t)/A  into  Equation  (2-27)  yields  the  PDE 


dC 

dt 


Q(t)  5C 


k C + D 


e2c 


A dz  dz 

The  design  of  feedback  controllers  for  PDEs  of  the  form  (2-28)  is  particularly 

challenging.  It  is  therefore  desirable  to  approximate  (2-28)  by  a set  of  ordinary 

differential  equations  (ODEs).  One  possible  approach  is  to  divide  the  volume  of  the 


(2-28) 
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contact  area  into  N smaller  volumes,  AVca  — Vca/N,  as  shown  in  Figure  2-2,  where  each 
volume  is  considered  a chemical  reactor. 


q,  c,  c,  c,  c,+1  cL 


Q(t) 

1 1 1 1 

Q(t) 

u(t) 

Dosage  Area 

1 1 Contact  Area  | | 



m 

Vda 

W 1 AVca  1 

II  II 

Z 


A Z1  ^ 3 3+1 

U L 

I 1 

Az 

Figure  2-2.  Schematic  of  the  wastewater  treatment  plant  when 
divided  into  N small  volumes  AVca  = Vca/N. 

th 

The  notation  Cj(t)  = C(t,zj)  represents  the  chlorine  concentration  of  the  i volume 

th 

section,  z\  is  the  horizontal  coordinate  that  denotes  the  end  of  the  i volume-section,  Az  = 

th 

Zi-Zi-1  is  the  length  of  the  i reactor,  and  AVca  = Vca  / N is  the  volume  of  each  reactor. 

All  other  quantities  are  as  given  in  Figure  2-1. 

Assuming  that  N is  large,  the  difference  equations 


dCj(t)  C.W-C.^t) 


i = 1,2,....N 


(2-29) 


dz  zj-zj_j 

can  be  used  to  approximate  the  spatial  partial  derivatives.  Combining  Equations  (2-28) 
and  (2-29)  yields  the  set  of  ordinary  differential  equations 


dCj(t) 

dt 


Q(t) 


A(zi  Zj,! ) 


(Ci(t)-Ci.1(t))- k^HC^t),  i = 1,2,..., N (2-30) 


Recognizing  that 
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A(zrzi_i)  = AVca 

produces  the  system  of  equations 


dCt(t) 

dt 

dC2(t) 

dt 


/ Q(t) 

AVca 

(Q(t) 

AVca 


)Ci(t)  + ( 
)C2(t)  + ( 


Q(t) 

AVca 

Q(t) 

AVca 


)C0(t)-k 

OQW-k 


rxn 


rxn 


c i(t) 

C2(t) 


(2-31) 


(2-32) 


(2-33) 


dCN (0  = ~(^-)CN  (t)  + (-^-)CN.!  (t)  - krxnCN(t) 


dt 


AV„ 


AV„ 


(2-34) 


ca  ’ ca 

which  can  be  written  in  a more  compact  notation  as 


dCj(t) 

dt 


Q(t) 

V AVca 


+ k 


rxn 


ci  (0  + (“~~")ci-l  (0>  i = 1,2,...,N  (2-35) 


AV, 


ca 


Finally,  the  outlet  chlorine  concentration  in  the  contact  area  is  found  as 


CL(t)  = CN(t)  (2-36) 

Equation  (2-35)  can  interpreted  as  a cascade  of  continuously-stirred  tank  reactors 

(CSTRs)  where  C0  is  the  input,  Q(t)  is  a disturbance,  Q are  the  states  for  i = 1,2,...,N , 

and  Cp  is  the  output.  Although  simplification  is  gained  because  the  system  is  now 
composed  of  N ODEs  instead  of  one  PDE,  it  is  important  to  note  that  each  equation  is 
nonlinear  because  of  the  presence  of  bilinear  terms  of  the  form  Q(t)Cj(t)  and  Q(t)Cj_](t). 

Finally,  just  how  large  N needs  to  be  in  order  to  accurately  describe  a wastewater 
treatment  plant  is  investigated  in  Section  2.3.  Note  that  the  behavior  described  by  the 
dispersion  term  of  (2-28)  is  approximated  when  N is  a finite  number.  Therefore,  it  is 
ignored  in  (2-35). 
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2.2.5  Linearized  State-Space  Model  of  a Plug  Flow  Reactor 

With  the  PDE  (2-27)  cast  as  a series  of  ODEs  (2-35),  the  final  step  in  creating  a 
linear  state-space  description  of  the  contact  basin  is  linearization  about  a steady  state.  A 
standard  linearization  process  applied  to  (2-35)  yields 


dCj(t) 

dt 


( Q 

vAVca 


+ k 


C. 


i-l 


AV  AV 


Q(t), 


ca  / 


(2-37) 


i = 1,2,...,N 

where  the  deviation  variables 


C,  —Ci-Cj, 

i = 0,1,2,..., N 

(2-38) 

Q Q(t)-Q 

(2-39) 

are  defined  as  a function  of  the  steady  state  values  C0 , C, , C2 , . . . CN , and  Q that  make  the 

right-hand  side  of  (2-37)  identically  zero.  The  model  (2-37)  can  now  be  cast  into  the  Nx2 
linear  state-space  form  (Kailath,  1980) 


x = Ax  + Bd 
y = Cx  + Dd 

where 


(2-40) 

(2-41) 


and 


[Cl  C2  ...  CN] 

(2-42) 

<iT=[Co  Q] 

(2-43) 
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where  x esJ?N,y  e^,d  es^2,A  e9?NxN,B  e'iRNx2,C  es}]lxN,  and  D eiH2x2. 

2.3  Plant  Data  Fitting  and  Parameter  Estimation 

In  order  to  have  confidence  that  a simulation  including  simplifying  assumptions  is 
an  accurate  approximation  of  a real  plant,  it  is  appropriate  to  validate  it  with  real  plant 
data.  Based  upon  the  model  developed  in  Section  2.2,  two  parameters  remain  to  be 
investigated  thoroughly  against  experimental  data:  namely,  the  reaction  rate  constant, 

krxrf  ar*d  the  number  N of  CSTRs  needed  to  adequately  match  the  experimental  data. 

Since  only  two  parameters  are  unknown,  they  can  be  easily  identified  through  a 
computationally  inexpensive  exhaustive  numerical  search.  Consider  the  cost  function 


1 ^ 


JN(krxn)  = -Z[cL(ti;N,krxn)-c[xp(ti)r 


i=l 


(2-48) 
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where  Cijxp(ti ) represents  the  experimental  concentration  values  measured  at  the  plant 
outlet  point  z = L at  time  tj  for  i = l,2,...,n,  and  is  the  output 

concentration  predicted  by  the  model  (2-35)  if  the  nonlinear  model  is  used  or  by  (2-37)  if 
the  linear  model  is  used.  Note  that  CL(ti;N,krxn)  is  a function  of  k„n— the  greater  k„n, 

the  lower  the  resulting  outlet  concentration  for  any  arbitrarily  chosen  value  of  N.  The 
functional  (2-48)  can  be  minimized  to  find  the  optimum  value  of  k^  for  a given  model 
of  order  N,  i.e. 


k«n(N)  = argmin  JN(kran)  (2-49) 

Figure  2-3  plots  (2-48)  as  a function  of  candidate  values  of  k^,  with  each  curve 

representing  a different  value  of  N.  The  figure  shows  that  for  N > 50,  the  value 

kmi  = 1 09x1 0-4  s'1  is  a component  of  the  optimizing  solution  of  (2-49) . Furthermore, 
close  inspection  reveals  that  the  minimum  over  N is  realized  at  roughly  N = 1 00 . In 
other  words,  as  N is  increased  from  to  25  toward  1 00,  becomes  increasingly  smaller, 
but  then  it  begins  to  increase  after  N becomes  greater  than  1 00.  This  is  in  reasonable 
agreement  with  tracer  tests  performed  at  the  facility  (Barrios  et  al.,  2002). 

It  should  be  noted  that  it  is  assumed  throughout  the  simulation  study  that  these 
values  are  constants.  In  practice  they  are  disturbances,  as  they  are  not  controlled,  and 
may  or  may  not  change  value.  The  additional  complexity  introduced  into  the  model  by 
allowing  them  to  vary  is  not  necessary  to  show  that  traditional  control  configurations  are 
inadequate  to  control  this  type  of  system. 
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Figure  2-3.  Map  of  the  cost  function  for  candidate  values  of  k^n, 
with  each  line  representing  a different  number  of  reactors  N.  A 
minimum  is  realized  forN=100  and  k^flOO)  = 1 .09xl0~4  s'1. 

Figure  2-4  shows  the  result  of  a simulation  using  a portion  of  experimental  plant 
data  gathered  at  the  Kanapaha  Wastewater  Treatment  Plant  (Demir,  2001).  The  bottom 
graph  shows  the  value  of  the  disturbance  Q(t)  versus  time.  Initially,  the  flow  is  held 

nearly  equal  to  0.5  m /s,  but  is  suddenly  decreased  to  slightly  more  than  0.4  m /s  after 

approximately  5.25  hours.  The  data  is  noise  free  because  the  available  flowrate  values 
are  simply  the  setpoint  of  a flowrate  controller,  as  the  flowrate  is  not  actually  measured 
on-line.  The  top  graph  plots  three  chlorine  concentrations  versus  time.  The  first  plot, 

depicted  by  a dashed  line,  is  the  concentration  Cq(t)  measured  at  the  beginning  of  the 
contact  basin.  The  experiment  involves  intentionally  increasing  Cq  at  approximately  t = 
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5 hours,  and  then  decreasing  its  original  value  at  approximately  t = 5.7  hours.  The 
second  concentration  plotted,  shown  as  a dashed-dot  line,  is  the  measured  value  of 
chlorine  at  the  outlet,  C^xp(t).  Finally,  the  third  line  is  solid,  and  corresponds  to  the 

output  concentration  CL(t,N)  predicted  by  passing  the  values  of  Co(t)  and  Q(t)  through 
the  nonlinear  model  (2-35)-(2-36)  for  the  special  case  of  N = 100  and 
kS(10°)=1-09xl0‘4  s'1. 
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Figure  2-4.  Fit  of  plant  data  using  value  of  that  minimizes  the 
Equation  (2-49)  for  N=100;  it  is  1^(100)=  1.09x1  O'4  s'1. 

Figure  2-4  suggests  that  the  proposed  model  parameters  N = 100  and  k^(100)  = 

1.09x10  4 s 1 provide  an  acceptable  fit.  However,  a few  characteristics  may  require 
further  explanation.  If  a linear  model  is  used,  its  predictions  are  most  accurate  near  the 
steady  state  about  which  the  nonlinear  model  has  been  linearized.  The  further  it  strays 
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from  that  point,  the  larger  the  errors.  Therefore,  at  the  expense  of  additional  computation 
time,  the  nonlinear  model  was  used  for  fitting  the  data.  Also,  the  model  output  appears 
much  smoother  than  the  plant  input;  this  is  because  the  series  of  CSTRs  acts  as  a low- 
pass  filter,  and  because  the  model  does  not  include  the  measurement  noise  associated 
with  the  observation  of  Cp(t). 

In  fact,  the  real  plant  output  (dashed-dot  line  in  Figure  2-4)  seems  to  be  a delayed, 
scaled,  and  smoothed  copy  of  the  true  plant  input  (dashed  line  in  Figure  2-4),  indicating 
that  the  plant  does  not  exhibit  truly  ideal  PFR  behavior.  Rather,  it  suffers  from  a small 
degree  of  mixing,  which  tends  to  better  mimic  the  output  of  a large  series  of  CSTRs.  In 
light  of  this,  instead  of  using  as  many  reactors  as  possible  to  minimize  in  an  attempt 
to  approximate  a PFR,  Figures  2-3  and  2-4  imply  that  there  is  a finite  number  of  reactors 
N that  best  describes  plant  data.  Figure  2-5  is  constructed  to  determine  the  optimal 
number  of  reactors  N that  minimizes  the  minima  of  the  curves  in  Figure  2-3.  It  plots  the 

number  of  reactors  N versus  JN  = (krxn)>  where  JN  is  given  in  (2-48),  k£J*  is  given  in 

(2-49),  and  focuses  on  a neighborhood  of  N = 100. 

It  is  clear  that  the  function  is  minimized  very  near  to  N - 100.  In  fact,  N = 102 
does  minimize  the  cost  function  to  a resolution  of  two  CSTRs.  As  with  any  data  fitting 
exercise,  higher  degrees  of  precision  can  be  found  with  ever  more  exhaustive  searches. 
However,  there  is  a limit  to  their  validity  beyond  a couple  of  decimal  places.  Therefore, 

the  values  N = 100  and  k^(100)  = 1.09x1  O'4  s'1  are  chosen  as  sufficient  for  the 
purposes  of  modeling  both  open-loop  and  closed-loop  performance. 
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Figure  2-5.  Optimization  of  N to  minimize  the  error  between  the 
nonlinear  model  and  the  true  plant  output. 


2.4  Control  Design  and  Assessment  of  Performance 

With  a viable  model  of  a wastewater  chlorination  system,  including  values  for  its 
parameters,  characteristics  of  the  particular  system  must  be  considered  in  selecting  an 
appropriate  control  scheme.  For  example,  the  system  is  affected  by  two  primary 
disturbances,  ammonia  chlorine  demand  and  flowrate.  Also,  there  is  a large  time  delay 
associated  with  the  disinfection  process,  implying  that  deadtime  compensation  of  some 
type  could  be  beneficial. 

Unlike  more  complicated  disturbances  such  as  flowrate  through  a fixed  volume 
(discussion  to  follow),  a cascade  control  system  is  a common  strategy  employed  to  reject 
load-type  input  disturbances  such  as  ammonia  chlorine  demand.  Therefore,  a cascade 
control  scheme  is  considered  in  addition  to  standard  feedback  design.  As  stated 
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previously,  the  dosing  area  portion  of  the  reactor  does  not  suffer  from  large  time  delay;  it 
is  the  contact  area  that  has  a large  deadtime.  Therefore,  a proportional-integral  (PI) 
controller  should  be  sufficient  to  achieve  satisfactory  control  in  the  slave  loop,  when  it  is 
utilized.  However,  the  extremely  large  delay  of  the  contact  basin  suggests  that  a 
deadtime  compensation  strategy,  in  addition  to  a standard  PI  controller,  should  be 
considered  for  the  master  controller.  (Note:  while  the  lone  controller  in  a standard 
feedback  system  is  generally  referred  to  as  a “feedback”  controller,  it  is  also  referred  to  as 
a “master”  controller  in  this  study  for  brevity.)  A Smith  Predictor  is  chosen  as  the 
deadtime  compensator.  Thus,  there  are  three  variables  pertinent  to  control  design,  and 
each  has  two  possible  values;  they  are  summarized  in  Table  2-3. 


Table  2-3.  Summary  of  the  total 
variations  to  consider  in  choosing  a 
control  design. 


| Control  Scheme 

1.  Standard  Feedback 

2.  Cascade  Feedback 

| Master  Controller 

1.  PI 

2.  PI  with  a Smith  Predictor 

1 Disturbances 

1 . Ammonia  Chlorine  Demand 

2.  Flowrate 

Because  the  effect  of  the  ammonia  chlorine  demand  on  the  system  is  fully 
realized  in  the  dosing  area  (fast  reactions  consume  chlorine  before  leaving  the  dosing 
area),  it  follows  that  a secondary  loop  (cascade  control)  should  be  effective  in  the 
rejection  of  this  disturbance.  Consequently,  cascade  control  is  anticipated  to  be  more 
effective  than  standard  feedback  control,  regardless  of  which  master  controller  is  chosen; 
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results  verifying  this  assertion  are  given  in  Appendix  B.  The  block  diagram  and  process 
used  for  obtaining  open-loop  responses,  controller  tunings,  and  simulation  results  in 
standard  feedback  configuration  is  analogous  to  that  discussed  in  detail  in  Section  2.4.2 
for  a cascade  configuration.  Therefore,  cascade  control  is  assumed  to  solve  the  problem 
of  chlorine  demand  as  a control  design  issue. 

In  addition,  flowrate  can  easily  be  measured,  and  input  concentration  is  often 
maintained  with  a (feedforward)  ratio  controller,  as  is  always  the  case  in  this  study. 

When  such  a strategy  is  employed,  the  ratio  controller  can  be  inserted  into  the  control 
loop  independent  of  all  other  design  choices.  For  the  purposes  of  modeling,  the  valve 
dynamics  are  assumed  to  be  much  faster  than  the  transport  delays  associated  with  a 
reactor  of  this  size.  As  a result,  the  model  of  a ratio  controller  is  omitted  within  the 
simulations  for  simplicity,  and  it  is  assumed  that  the  desired  input  concentration  u(t)  can 
be  achieved  instantaneously  (i.e.,  the  ratio  controller  has  no  dynamics;  rather,  it  is  simply 
a ratio  station  that  performs  a calculation).  In  practice,  any  error  introduced  by  using  a 
feedforward  controller  is  automatically  corrected  by  the  slave  controller  that  is  being  used 
to  reject  the  ammonia  chlorine  demand  disturbance,  thus  reinforcing  the  argument  for 
using  cascade  control  rather  than  standard  feedback  control. 

However,  input  concentration  is  only  one  of  two  ways  in  which  flowrate  effects 
the  system;  residence  time  is  the  other.  Residence  time  is  an  effect  seen  throughout  the 
reactor,  so  it  is  a characteristic  that  cannot  be  compensated  for  solely  by  a secondary 
(slave)  loop.  The  control  design  problem  is  now  to  select  the  proper  master  controller  (to 
be  used  in  a cascade  control  scheme)  based  upon  the  closed-loop  performance  of  the 
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system  when  subject  to  realistic  flowrate  patterns.  In  other  words,  the  pertinent  control 
design  issues  presented  in  Table  2-3  are  reduced  to  those  given  in  Table  2-4. 


Table  2-4.  Summary  of  the  simplified 
total  variations  to  consider  in  choosing  a 
control  design. 


Control  Scheme 

Cascade  Feedback  • 

Master  Controller 

1.  PI 

2.  PI  with  a Smith  Predictor 

1 Disturbances 

Flowrate 

A block  diagram  of  the  cascade  control  scheme  is  shown  in  Figure  2-6,  where 
Gsc,  Gmc,  and  RC  are  the  slave,  master,  and  ratio  controllers,  respectively.  In  addition, 

Gda  is  the  dosage-area  model,  Gca  is  the  contact-area  model,  rsc  and  rmc  are  the  slave  and 


master  setpoints,  respectively,  and  esc  and  emc  are  the  slave  and  master  error  signals, 

respectively.  Finally,  Rd(t)  = Qd(0/Q(0  is  the  flowrate  ratio  required  of  the  ratio 
controller.  All  other  quantities  are  the  same  as  defined  in  Section  2.2. 


w 

Q(t) 

5(t) 


Figure  2-6.  Block  diagram  of  a wastewater  treatment  plant 
under  cascade  control. 
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An  input  of  chlorine  dosing  stream  flowrate  Qd(t)  (and  with  constant  chlorine 

concentration  Cj)  is  received  in  the  dosing  area  from  the  ratio  controller  RC.  The 

chlorine  is  mixed  with  the  incoming  untreated  water  arriving  at  a flowrate  Q(t),  and 
chlorine  consumption  due  to  reaction  with  ammonia  occurs  to  a significant  extent.  Then, 
the  chlorinated  stream  passes  from  the  dosing  area  into  the  contact  area  transported  at  the 

same  flowrate  Q(t)  and  with  an  established  initial  concentration  Co(t). 

It  is  during  the  extended  period  of  time  that  the  liquid  is  in  the  contact  area  that 
the  disinfecting  and  slow  reactions  take  place,  as  well  as  volatilization  and  UV  reactions. 

The  inlet  concentration  Co(t)  is  reduced  to  the  outlet  concentration  value  Cp(t)  during  the 
time  the  flow  traverses  the  length  of  the  reactor.  The  exit  concentration  Cp(t)  is 
measured  and  compared  to  the  master  setpoint  rmc.  This  signal  is  processed  in  the  master 
controller  Gmc  which  prescribes  an  appropriate  slave  setpoint  rsc.  The  latter  is  compared 
to  the  measurement  Co(t)  to  produce  the  feedback  error  esc,  which  in  turn  is  processed  by 

the  slave  controller  Gsc  producing  the  prescription  Rd(t),  which  is  sent  to  the  ratio 
controller.  The  ratio  controller  uses  the  information,  along  with  a measurement  of  Q(t), 
to  produce  Qd(t),  thereby  completing  the  loop. 

2.4.1  Description  and  Tuning  of  the  Slave  PI  Controller 

The  slave  proportional-integral  controller  (Ogunnaike  and  Ray,  1994)  is  of  the 


standard  form 
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Gsc  (s)  = Ksc  ( 1 -t  — ) (2-50) 

Tscs 

where  Ksc  is  the  gain  of  the  controller,  and  isc  is  the  integral  time  constant  of  the 
controller. 

2.4. 1.1  Open-loop  response  of  the  dosing  area  model  to  step  changes  in  dosing  rate 
and  tuning  parameters  of  a slave  PI  controller 

In  cascade  systems,  the  inner  loop  (slave  loop)  is  normally  more  aggressive,  and 
must  be  tuned  first  (Ogunnaike  and  Ray,  1994).  The  outer  loop  is  left  open  during  this 
procedure  (i.e.  the  gain  of  the  master  controller  is  set  equal  to  zero).  The  range  of  typical 

3 3 

(as  opposed  to  extreme)  flowrates  Q(t)  for  the  Kanapaha  facility  is  0.2  m /s  to  0.6  m /s 


0 0.05  0.1  0.15  0.2  0.25  0.3  0.35  0.4  0.45  0.5 


t [hrs] 

Figure  2-7.  Open-loop  response  of  the  dosing  area  for  a step  in 
dosage  u(t)  = Qd(t)Ccj  from  2.451  to  2.651  g/s  at  disturbance  value 
of  Q(t)  = 0.2  m3/s. 
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3 3 

(recall  that  the  nominal  value  is  0.4  m /s).  Accordingly,  0.6,  0.4,  and  0.2  m /s  are 

designated  as  high,  middle,  and  low  values  of  the  disturbance  Q(t).  Figure  2-7  shows  the 
open-loop  response  of  the  dosing  area  to  a step  change  in  input  concentration  u(t)  from 
12.25  ppm  to  13.25  ppm  at  time  0.1  hr  as  predicted  by  the  model  (2-15).  The  top  plot 

shows  the  output  concentration  Cg(t)  as  a function  of  time,  the  center  plot  shows  the  input 

concentration  u(t)  as  a function  of  time,  and  the  bottom  plot  shows  that  the  flowrate  is 

3 

held  constant  at  the  value  Q = 0.2  m /s  (the  ammonia  chlorine  demand  8 is  also  held 
constant  at  3 ppm,  but  not  depicted). 

The  first-order  plus  deadtime  (FODT)  fit  (Seborg,  Edgar,  and  Mellichamp,  1989) 


Figure  2-8.  FODT  fit  of  the  dosing  area  for  a step  in 
Qd(t)Cd  from  2.451  to  2.651  g/s;  the  observed  parameters 
are  Kda  = 5.0  ppm-s/g,  xda  = 280  s,  and  0da  = 0 s. 
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ol  the  concentration  output  Co(t)  of  the  dosing  area  for  the  open-loop  step  change  in 

dosage  input  u(t)  from  2.451  g/s  to  2.651  g/s  (the  top  plot  in  Figure  2-7)  is  shown  in 
Figure  2-8.  Output  from  the  simulation  is  represented  with  a solid  line,  and  the  FODT 
model  is  represented  with  a dashed-dot  line.  The  criterion  used  to  decide  which  curve 
yields  the  best  fit  is  minimization  of  the  sum  of  the  squared  error  (SSE). 

The  procedure  is  repeated  with  a step  down  in  input,  and  at  middle  and  high  flow 
conditions,  resulting  in  a total  of  six  sets  of  FODT  parameters.  They  are  then  averaged, 
and  the  averages  are  used  to  tune  the  slave  controller.  Table  2-5  lists  the  inputs  and 
outputs  for  each  trial,  the  observed  FODT  transfer  function  parameters  for  each  case,  and 
the  averages  of  the  values.  The  first  two  columns  report  the  value  of  the  input  u(t)  in  g/s 
before  and  after  the  step  change,  respectively.  The  next  two  columns  give  the  value  of 

the  output  Co(t)  in  ppm  both  before  the  step  and  after  the  system  has  been  allowed  to 


Table  2-5.  Dosing  area  transfer  function  parameters  for 
various  steps  in  control  variable  u(t)  = Qd(t)Cd  for 
various  values  of  disturbance  Q(t). 


Initial 

Final 

Initial 

Final 

Plant 

Kria 

Tda 

^da 

S.S.  DA 

S.S.  DA 

S.S.  DA 

S.S.  DA 

Flowrate 

Input 

Input 

Output 

Output 

Q(t) 

Qd(t)Cd 

Qd(t)Cd 

c0(t) 

c0(t) 

[s] 

[S[ 

[g/sl 

[g/s] 

[ppm] 

[ppm] 

[m  /s] 

[ppm-s/g] 

2.451 

2.651 

9.25 

10.25 

0.2 

5.0 

280 

0 

2.451 

2.251 

9.25 

8.26 

0.2 

5.0 

280 

0 

2.926 

3.326 

4.31 

5.31 

0.4 

2.5 

140 

0 

2.926 

2.526 

4.31 

3.31 

0.4 

2.5 

140 

0 

1 3.805  i 

4.405 

3.34 

4.34 

0.6 

1.7 

100 

0 

3.805 

3.205 

3.34 

2.34 

0.6 

1.7 

100 

0 

Average 

3.1 

170 

0 
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reach  its  new  steady-state,  respectively.  The  fifth  column  gives  the  value  of  the  flowrate 
disturbance  Q(t)  for  each  particular  trial,  and  the  final  three  columns  report  the  observed 

FODT  parameters,  Kda,  ?da>  and  9da.  Finally,  average  vales  for  each  of  the  parameters 

are  listed  across  the  bottom  row. 


Table  2-6  repeats  the  average  values,  and  reports  the  resulting  values  of  slave 
controller  gain  Ksc  and  integral  time  constant  xsc.  In  this  case,  they  are  found  by  the 

trial-and-error  method,  since  the  ratio  of  0 / x does  not  fall  into  the  suggested  range  of 
0. 1 < 0 / x < 1 .0  set  forth  by  standard  tuning  rules  such  the  ITAE  and  ISE  methods 
(Ogunnaike  and  Ray,  1994),  primarily  because  the  FODT  fit  of  an  ideal  CSTR  model 
results  in  a zero  delay  value. 


Table  2-6.  Estimated  dosing  area  transfer  function 
parameters  and  resulting  slave  controller  parameters 
found  by  trial-and-error  tuning. 


Transfer  Function 

Numerical  Value 

Controller 

Numerical  Value 

i Parameters 

Parameters 

Kda 

3.1  ppm-s/g 

Ksc 

0.025  g/ppm-s 

* 

Tda 

170s 

"frc 

50  s 

9da 

Os 

2.4. 1.2  Closed-loop  response  of  the  dosing  area  model  for  step  changes  in 
concentration  setpoint 

Figure  2-9  shows  three  closed-loop  responses  of  the  dosing  area  model  (2-15)  for 
a setpoint  change  in  rsc  from  2 to  3 ppm  at  time  0.1  hr,  which  is  depicted  as  a dotted  line. 

The  top  plot  shows  the  output  concentration  Cg(t)  as  a function  of  time,  the  center  plot  is 


the  input  to  the  dosing  area  u(t)  administered  by  the  slave  controller  as  a function  of  time. 
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and  the  bottom  plot  is  the  water  flowrate  Q(t)  through  the  dosing  area  as  a function  of 

_ 3 

time.  The  dashed-dot  line  represents  the  result  for  a low-flow  disturbance  (0.2  m /s),  the 

3 

solid  line  represents  the  calculation  for  a medium-flow  disturbance  (0.4  m /s),  and  the 

3 

dashed  line  represents  the  response  for  a high-flow  disturbance  value  (0.6  m /s). 

The  controller  accomplishes  a setpoint  tracking  change  for  each  disturbance  value 
almost  instantly;  in  addition,  overshoot  is  minimal,  and  oscillations  are  almost  non- 
existent. The  primary  reason  is  that  the  dosing  area  is  based  upon  an  ideal  CSTR,  which 
has  no  time  delay  from  the  input  to  the  output.  In  reality,  a dosing  area  is  likely  to  have 
some  small  delay.  However,  note  that  the  slave  controller  accomplishes  the  setpoint 
change  within  a very  small  fraction  of  an  hour.  Even  with  a small  delay  and  the 
associated  dynamics,  the  closed-loop  response  would  resemble  a step  when  viewed  on 
the  time  scale  of  the  contact  area.  Therefore,  the  model  provides  an  acceptable 
representation  of  a slave  controller  implemented  around  a real  dosing  area. 

In  addition  to  setpoint  tracking,  the  slave  loop  is  required  to  reject  disturbances  in 
ammonia  chlorine  demand.  The  simulation  studies  show  that  here  again,  the  performance 
is  comparatively  fast,  and  the  closed-loop  response  essentially  appears  to  be  a flat  line 
when  viewed  on  the  time  scale  of  the  contact  area.  These  responses  are  not  shown  for 
brevity.  Note  also  that  the  simulations  are  performed  with  the  slow  reactions  included, 
yet  the  steady-state  behavior  of  the  dosing  area  is  virtually  identical  at  all  flowrates, 
thereby  confirming  the  assumption  that  only  the  fast  reactions  are  important  in  this  part  of 


the  reactor. 
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Figure  2-9.  Closed-loop  response  of  the  dosing  area  using  a PI 

controller  with  Ksc  = 0.025  g/ppm-s  and  rsc  =50  s for  disturbance 

3 3 3 

values  Q(t)  = 0.2  m /s,  0.4  m /s,  and  0.6  m Is. 


2.4.2  Description  and  Tuning  of  the  Master  PI  Controller 

With  the  inner  loop  tuned,  the  outer  master  loop  must  be  tuned  next.  A PI 

controller  is  selected  for  Gmc  of  the  familiar  form 


Gmc(s)  = Kmc(l+— ) 

Tmcs 

where  Kmc  is  the  controller  gain  and  Tmc  is  the  controller  integral  time  constant. 


(2-51) 
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2.4.2.1  Open-loop  response  of  the  contact  area  model  to  step  changes  in  input 
concentration  and  tuning  parameters  of  a master  PI  controller 

With  the  inner  slave  loop  operating,  a step  change  in  the  variable  C0(t)  is  passed 

through  the  entire  system  to  obtain  an  open-loop  response  for  the  contact  area.  Figure  2- 

10  shows  the  results.  The  top  plot  shows  the  output  concentration  CL(t)  as  a fiinction  of 

time,  the  center  plot  shows  that  the  setpoint  C0(t)  (i.e.  rsc)  is  changed  from  9.25  ppm  to 
10.25  ppm  at  time  2 hr,  and  the  bottom  plot  shows  that  the  disturbance  Q(t)  is  held 

3 

constant  at  Q = 0.2  m /s.  In  addition,  ammonia  chlorine  demand  is  again  held  constant  at 
5 = 3 ppm,  but  not  shown. 


t [hr] 


Figure  2-10.  Open-loop  response  of  the  contact  area  with  the  slave 
controller  in  operation  about  the  dosing  area. 

Inspection  of  the  first  two  graphs  reveals  that  setpoint  adjustments  of  the  slave 
controller  cause  changes  in  the  output  variable  CL(t)  that  are  representative  of  a very  high 
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order  system.  However,  the  response  can  be  effectively  approximated  with  a FODT 
model.  The  results  of  the  FODT  fit  are  in  turn  used  to  tune  the  master  controller. 

The  FODT  fit  of  the  concentration  output  Cp(t)  of  the  contact  area  for  the  open- 

loop  step  change  in  input  C0(t)  from  9.25  ppm  to  10.25  ppm  (the  top  plot  in  Figure  2-10) 

is  shown  in  Figure  2-11.  Output  from  the  simulation  is  represented  with  a solid  line,  and 
the  FODT  model  is  represented  with  a dashed-dot  line.  The  criterion  used  to  decide 
which  curve  gives  the  best  fit  is  again  minimization  of  the  SSE. 


Figure  2-11.  FODT  fit  of  the  contact  area  (while  the  slave 
controller  is  in  operation  around  the  dosing  area)  for  step  change  in 

Co(t)  from  9.25  ppm  to  10.25  ppm  with  the  flow  held  constant  at  Q 
= 0.2  m /s;  the  observed  parameters  are  Kca  = 0.216  ppm-s/g,  rca 
= 1650  s,  and  0ca  = 12510  s. 

When  the  process  of  obtaining  an  open-loop  response  and  a FODT  model  fit  is 
repeated  for  several  operating  conditions,  the  resulting  values  are  compared  in  an  effort  to 
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find  parameters  most  representative  of  the  system  at  all  conditions.  The  values  found  to 
minimize  the  error  for  various  operating  conditions  are  listed  in  Table  2-7.  The  first  two 
columns  represent  the  beginning  and  ending  steady-state  inputs,  and  the  next  two 
columns  report  the  beginning  and  ending  steady-state  outputs.  The  fifth  column  gives  the 
plant  disturbance  Q(t)  value  during  the  trial.  Finally,  the  last  three  columns  list  the  values 

of  Kca,  Tea,  and  9ca  that  minimize  the  error  of  the  FODT  model,  respectively.  The 

averages  of  all  five  trials  are  found  across  the  bottom  row.  (Note  that  since  steps  up  and 
steps  down  yield  the  same  results  for  a given  disturbance,  more  flowrate  operating  points 
are  considered  in  lieu  of  redundant  information.) 


Table  2-7.  Contact  area  transfer  function  parameters  for 

various  steps  in  input  variable  Co(t)  for  various  values  of 
disturbance  Q(t). 


Initial 

S.S. 

Input 

c0(t) 

[ppm] 

Final 

S.S. 

Input 

c0(t) 

[ppm] 

Initial 

S.S. 

Output 

c0(t) 

[PPm] 

| Final 
S.S. 
Output 

c0(t) 

[PPm] 

Plant 

Flowrate 

Q(t) 

[m3/s] 

Kca 

fPPm-s/g] 

Tca 

Is] 

®ca 

[s] 

9.25 

10.25 

2.00 

2.22 

0.2 

0.216 

1650 

12510 

5.57 

6.57 

2.00 

2.36 

0.3 

0.359 

1100 

8390 

4.31 

5.31 

2.00 

2.46 

0.4 

0.464 

820 

6310 

3.70 

4.70 

2.00 

2.54 

0.5 

0.540 

660 

5060 

3.34 

4.34 

2.00 

2.60 

0.6 

0.599 

550 

4220 

Average 

0.436 

960 

7300 

Using  the  average  values  for  the  contact  area,  it  follows  that  the  ratio  0 / x again 


does  not  fall  into  the  standard  range  0. 1 < 0 / x < 1 .0.  Asa  result,  the  controller  is  tuned 
as  before  via  trial-and-error.  Finally,  recall  that  because  the  master  controller  must 
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primarily  maintain  a value  of  output  concentration  Cp(t),  it  is  tuned  for  optimal 
regulatory  performance,  rather  than  setpoint  tracking. 


Table  2-8.  Estimated  contact  area  transfer  function 
parameters  and  resulting  master  controller  disturbance- 
rejection  parameters  found  by  trial-and-error  tuning. 


Transfer  Function 
Parameters 

Numerical  Value 

Controller 

Parameters 

Numerical  Value 

0.436  ppm-s/g 

Kmc 

0.6  g/ppm-s 

* 

Tca 

960  s 

fine 

2500  s 

<4 

7300  s 

The  combined  performance  of  the  two  controllers  is  assessed  next  when  applied 
to  the  linear  plant  model  for  various  disturbance  Q(t)  conditions.  Even  though  the 


controller  designs  must  ultimately  be  applied  to  the  more  realistic,  nonlinear  plant  model, 
the  linear  closed-loop  response  is  important  for  several  reasons.  First,  the  controller  that 
is  being  applied  to  the  system  is  linear;  hence,  if  it  cannot  yield  satisfactory  results  for  the 
linear  plant  model,  then  it  is  extremely  unlikely  that  it  would  be  effective  on  the  original 
nonlinear  model.  In  addition,  because  software  used  for  simulation  (Matlab)  is  matrix- 
based,  it  can  simulate  a state-space  model  ten  to  one  hundred  times  faster  than  a user- 

defined  function.  This  is  extremely  useful  in  areas  of  debugging,  tuning,  and  preliminary 
validation. 

2.4.2.2  Closed-loop  response  of  the  entire  model  for  step  changes  in  flowrate 

After  the  tuning  choices  are  tested  and  verified  on  the  linear  model,  they  are  used 
to  attempt  to  control  the  nonlinear  model.  Figure  2-12  records  the  response  of  the  system 
when  subjected  to  a small  early  step-up  in  disturbance  Q(t)  at  time  2 hr,  followed  by  a 
large  step-down  in  disturbance  at  time  12  hr,  and  finally  another  step-up  at  time  34  hr, 
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returning  to  the  original  flowrate  Q(t).  The  first  plot  shows  the  output  concentration  for 
the  linear  plant  (dashed-dot  line),  the  nonlinear  plant  (solid  line),  and  the  setpoint  (dotted 
line)  as  functions  of  time.  Because  regulations  require  a minimum  of  1 ppm  chlorine  in 
water  discharged  from  the  plant  (Florida  Department  of  Environmental  Protection,  1996), 
a value  of  2 ppm  is  chosen  for  the  setpoint,  providing  a margin  of  safety.  The  second  plot 

shows  the  chlorine  concentration  entering  the  contact  area  C0(t)  for  each  of  the  two 

models  over  time.  Finally,  the  disturbance  Q(t)  is  shown  in  the  third  plot  as  a function  of 
time. 


t [hr] 

Figure  2-12.  Closed-Loop  response  of  linear  and  nonlinear 
models  with  a pair  of  controllers  arranged  in  cascade 
configuration  and  subjected  to  a step  changes  in  flowrate. 

The  well-tuned  controller  accomplishes  an  excellent  response  for  both  the  linear 
and  nonlinear  plant.  As  expected,  it  performs  better  on  the  linear  plant,  returning  the 
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output  concentration  Cp(t)  to  the  desired  value  within  eight  hours,  even  for  the  large  step 

change  in  disturbance  Q(t)  at  time  12  hr.  The  large  step  is  designed  to  represent  a 
condition  of  low  residence  time  (less  than  1 */2  hours)  immediately  followed  by  a 
condition  of  high  residence  time  (approximately  4 hours),  creating  a worst-case  scenario. 
As  a result,  eight  hours  (two  residence  times)  is  not  an  unacceptably  large  amount  of  time 
for  the  master  controller  to  bring  the  system  back  to  the  setpoint.  Furthermore,  careful 
selection  of  the  tuning  parameters  also  results  in  a nonlinear  closed-loop  response  with 
acceptable  performance,  almost  matching  that  of  the  closed-loop  linear  system. 

2. 4. 2. 3 Closed-loop  response  of  the  entire  model  for  a dynamic  flowrate 

Figure  2-13  shows  the  performance  of  the  same  control  systems  applied  to  the 
entire  plant  model  when  the  flowrate  is  modeled  as  sinusoidal  wave,  which  is  chosen  to 
approximate  a smooth,  cyclical  water  consumption  profile  of  a typical  day.  The  sine 
wa\e  begins  at  time  2 hr,  its  period  is  24  hr,  and  the  amplitude  is  three-fourths  of  the 

nominal  flowrate  (0.4  m /s),  or  0.3  m3/s.  The  result  is  a disturbance  Q(t)  that  varies  from 
3 3 

0.1  m /s  to  0.7  m /s.  These  values  are  chosen  because  they  represent  a 7:1  ratio  of 

maximum-to-minimum  flowrate,  which  is  not  an  unrealistic  occurrence. 

The  linear  plant  is  controlled  with  acceptable  performance,  experiencing  only 

low-amplitude  fluctuations  in  output  CL(t)  that  lag  predictably  behind  the  flowrate 

profile.  However,  it  is  clear  that  a sinusoidal  flow  of  significant  amplitude  severely 
degrades  the  quality  of  the  controller  performance  for  simulations  of  the  nonlinear  plant. 

Figure  2-1 3 shows  that  output  concentration  values  CL(t)  are  as  high  as  6.0  ppm  and  as 

low  as  0.5  ppm.  The  conclusion  is  that  cascade  control  of  a wastewater  treatment  plant 
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Figure  2-13.  Closed-Loop  response  of  linear  and  nonlinear 
models  with  a pair  of  controllers  arranged  in  cascade 
configuration  subjected  to  varying  flow. 

with  a PI  slave  controller  and  a PI  master  controller  is  viable  if  the  plant  is  linear,  but  it  is 
inadequate  in  more  realistic  circumstances  that  include  the  nonlinear  dynamics  of  the 
system. 

2.4.3  Description  and  Tuning  of  the  Master  PI  Controller  and  Smith  Predictor 

A PI  master  controller  with  a Smith  Predictor  (Ogunnaike  and  Ray,  1994)  can  be 
used  in  an  attempt  to  reduce  the  effect  of  the  delay  and  improve  performance.  The 
structure  of  a Smith  Predictor  applied  to  a PI  master  controller  is  given  by  the  expression 


K-mcO  + ) 


Gmc(s)  = 


Tmcs 


‘mcv  1 ;*^ca  — * 0~e  0caS) 

Tmcs  Tcas+1 


1 + KmcO+ )KC 


(2-52) 
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where  the  numerator  factors  Kmc(l  h — ? — ) represent  a standard  PI  controller,  K*  is  the 

-T  Q ~ td 

lmc5 

anticipated  value  of  the  contact  area  gain,  x*a  is  the  anticipated  value  of  the  contact  area 
time  constant,  and  0ca  is  the  expected  delay  of  the  contact  area. 

As  in  the  case  of  simple  PI  master  controller,  the  tuning  parameters  Kmc  and  rmc 

must  be  selected  to  achieve  optimum  performance.  With  the  effect  of  the  delay  reduced, 
however,  the  controller  tunings  change,  and  may  become  more  aggressive.  Furthermore, 
with  the  delay  theoretically  removed  from  the  model,  the  trial-and-error  tuning  method  is 
again  the  most  reliable.  Using  the  same  FODT  model  for  the  contact  area  found  in  Table 

2-7,  the  new  values  of  Kmc  and  xmc  are  reported  in  Table  2-9.  Note  that  K*a,  x*a,  and 
* 

0ca  are  simply  the  contact  area  model  parameters  found  in  Table  2-7. 


Table  2-9.  Contact  area  transfer  function  parameters  and 
resulting  disturbance  rejection  trial-and-error  tunings 
when  a Smith  Predictor  is  incorporated. 


Transfer  Function 
Parameters 

Numerical  Value 

Controller 

Parameters 

Numerical  Value 

Kc*a 

0.436  ppm-s/g 

^mc 

0.6  g/ppm-s 

* i 

Tca 

960  s 

Tmc 

1000  s 

Qca 

7300  s 

2.4.3.1  Closed-loop  response  of  the  entire  model  for  step  changes  in  flowrate  using  a 
Smith  Predictor 

Figure  2-14  shows  the  simulation  results  for  the  system  under  PI  control  with  a 
Smith  Predictor  subject  to  three  separate  constant  delays.  The  flowrate  profile  (three  step 
changes)  and  setpoint  (held  constant  at  2 ppm)  are  identical  to  that  used  in  Figure  2-12, 
using  the  same  line-type  convention  as  well.  It  is  anticipated  that  the  response  is 
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improved  over  the  case  of  a standard,  single  PI  controller,  provided  that  the  model  values 
* * * 

Kca>  Tca  and  0Ca  are  sufficiently  accurate. 


t [hr] 

Figure  2-14.  Closed-Loop  response  of  linear  and  nonlinear 
models  with  a pair  of  controllers  and  a Smith  Predictor 
arranged  in  cascade  configuration  subjected  to  a step 
changes  in  flowrate. 

For  these  conditions,  the  cascade  control  configuration  with  a Smith  Predictor 


added  to  the  PI  controller  in  the  master  loop  yields  an  almost  identical  response  to  that 

with  a standard  PI  in  the  master  loop,  both  for  the  linear  and  nonlinear  cases. 

2.4.3.2  Closed-loop  response  of  the  entire  model  for  a dynamic  flowrate  using  a 
Smith  Predictor 

When  the  flowrate  profile  Q(t)  of  the  system  under  cascade  PI  control  with  the 
addition  of  a Smith  Predictor  to  the  master  loop  is  changed  to  a sine  wave,  the  benefit  of 
the  Smith  Predictor  is  again  unapparent,  as  shown  in  Figure  2-15  Similarly  to  Figure  2- 
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1 1,  the  periodic  flowrate  Q(t)  profile  begins  at  time  2 hr,  has  a period  of  24  hr,  and  has  an 

3 

amplitude  of  0.3  m /s. 


Figure  2-15.  Closed-Loop  response  of  linear  and  nonlinear 
models  with  a pair  of  controllers  and  a Smith  Predictor 
arranged  in  cascade  configuration  subject  to  varying  flow. 

In  the  continuously  varying  flowrate  scenario,  only  the  linear  plant  simulation 
produces  a satisfactory  response.  As  in  Figure  2-11,  the  output  is  a smooth  sine  wave 
that  lags  behind  the  disturbance.  However,  the  controller  fails  almost  completely  to 
counter  the  effects  of  continuously  varying  delay  when  applied  to  the  nonlinear  plant,  and 
the  output  is  virtually  the  same  as  that  obtained  without  a Smith  Predictor  under  varying 
flow  conditions.  The  result  is  in  fact  slightly  worse  than  for  simple  cascade  control,  since 
the  given  that  the  output  ranges  from  0.5  ppm  to  6.4  ppm. 
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Because  the  Smith  Predictor  is  generally  considered  effective  only  within  the 
range  ±30%  of  the  estimated  delay  (Seborg,  Edgar,  and  Mellichamp,  1989),  and  given 
that  the  deadtime  can  often  vary  five-to-tenfold  in  a wastewater  treatment  plant,  the  lack 
of  performance  of  the  Smith  Predictor  observed  is  easily  explained.  Although  a standard 
PI  controller  is  adequate  for  the  slave  loop,  a PI  controller  either  with  or  without  a Smith 
Predictor  should  not  be  used  as  a master  controller.  A more  sophisticated  controller  form 
is  needed  to  achieve  satisfactory  performance. 

2.4.4  Description  and  Tuning  of  the  Master  PI  Controller  and  Adaptive  Smith 
Predictor 

The  response  of  the  system  under  a PI  master  controller  with  a Smith  Predictor 
(Ogunnaike  and  Ray,  1994)  is  not  adequate,  at  least  partly  because  the  value  of  the  delay 
0ca  changes  dramatically  from  one  operating  point  to  the  next.  This  can  be  seen  in  Table 
2-7,  which  lists  the  FODT  fit  parameters  of  the  model  for  open-loop  steps  changes  in 
input  concentration  C0  at  various  constant  disturbance  values  Q.  Note  also  that  the 

values  of  the  other  parameters,  Kca  and  xca , vary  with  operating  point.  As  a result,  it  is 

clear  that  a Smith  Predictor  with  parameter  estimates  at  the  nominal  flowrate  (i.e.  0.4 
3 

m /s  in  this  study)  will  experience  error  in  all  three  parameters  as  the  system  moves  to 

different  operating  points.  It  must  be  assumed  that  this  will  serve  to  only  exacerbate  the 
degradation  of  performance.  However,  each  parameter  follows  a trend,  and  the  three 
parameters  are  fit  extremely  well  by  the  equations 

K(Q)  = -1.407Q2  +2.073Q  — 0.140  (2-53) 


x(Q)  = (0.003Q- 0.000001  )•' 


(2-54) 
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0(Q)  = (0.0004Q-  0.00000 1)’1  (2-55) 

where  each  parameter  is  updated  continually  in  the  controller  as  a function  of  flowrate. 

This  vanation  of  a Smith  Predictor  is  an  adaptive  Smith  Predictor.  The  structure  of  a PI 

master  controller  with  an  adaptive  Smith  Predictor  is  given  by  the  expression 

K„(1+—  ) 

G~<s)  = j — ) (2-56) 

* + K,(l  + )K„x-!~7(l-e-1’-*) 

Tmcs  xC3s+l 

where  the  numerator  factors  Kmc(l  h — ? — ) represent  a standard  PI  controller,  K.  is  the 

Tmcs 

current  anticipated  value  of  the  contact  area  gain,  f ca  is  the  current  anticipated  value  of 
the  contact  area  time  constant,  and  0ca  is  the  current  expected  delay  of  the  contact  area. 

As  in  the  previous  cases,  the  tuning  parameters  Kmc  and  xmc  must  be  selected  to 

achieve  optimum  performance.  By  using  a dynamic  algorithm,  the  controller  tunings 
may  become  even  more  aggressive  than  for  the  case  with  a traditional  Smith  Predictor. 

The  FODT  model  for  the  contact  area  found  in  Table  2-7  is  again  used  as  the  starting 
point,  the  new  values  of  Kmc  and  xmc  are  reported  in  Table  2-10. 


Table  2-10.  Contact  area  transfer  function  parameters 
and  resulting  disturbance  rejection  trial-and-error  tunings 
when  a Smith  Predictor  is  incorporated. 


Transfer  Function 

Numerical  Value 

Controller 

Numerical  Value 

Parameters 

Parameters 

Kca 

0.436  ppm-s/g 

Kmc 

0.6  g/ppm-s 

Tca 

960  s 

Tmc 

200  s 

Qca 

7300  s 
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Note  that  controller  gain  does  not  increase,  as  is  the  case  for  many  deadtime-compensated 
systems,  but  the  integral  action  is  much  greater,  implying  that  there  is  potential  for 
improved  performance.  Interestingly,  these  parameters  are  valid  for  the  nonlinear  case 
only.  While  many  systems  tolerate  more  aggressive  tunings  for  the  linear  case,  the  PI 
master  controller  with  an  adaptive  Smith  Predictor  does  not.  In  order  to  keep  the  linear 
system  stable,  the  integral  time  constant  of  2000  s must  be  used. 

2.4.4. 1 Closed-loop  response  of  the  entire  model  for  step  changes  in  flowrate  using 

an  adaptive  Smith  Predictor 

Figure  2-16  shows  the  simulation  results  for  the  system  under  PI  control  with  an 
adaptive  Smith  Predictor  subject  to  the  same  three  separate  constant  delays  as  in  Section 


t [hr] 

Figure  2-16.  Closed-Loop  response  of  linear  and  nonlinear 
models  with  a pair  of  controllers  and  an  adaptive  Smith 
Predictor  arranged  in  cascade  configuration  subjected  to  a 
step  changes  in  flowrate. 
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2.4.3.  The  flowrate  profile  (three  step  changes)  and  setpoint  (held  constant  at  2 ppm)  are 
identical  to  that  used  in  Figure  2-12,  using  the  same  line-type  convention  as  well.  It  is 
anticipated  that  the  response  is  improved  over  the  cases  of  a standard,  single  PI 
controller,  as  well  as  that  of  the  PI  master  controller  with  a traditional  Smith  Predictor. 

For  these  conditions,  the  cascade  control  configuration  with  an  adaptive  Smith 
Predictor  added  to  the  PI  controller  in  the  master  loop  yields  a response  that  is 
significantly  better  than  with  a standard  PI  in  the  master  loop,  at  least  for  the  nonlinear 
case.  The  reason  is  likely  to  be  that  the  master  controller  tunings  for  simple  cascade  and 
cascade  with  a Smith  Predictor  cannot  tolerate  such  an  aggressive  integral  time. 

2.4.4.2  Closed-loop  response  of  the  entire  model  for  a dynamic  flowrate  using  a 

Smith  Predictor 

When  the  flowrate  profile  Q(t)  of  the  system  under  cascade  PI  control  with  the 
addition  of  an  adaptive  Smith  Predictor  to  the  master  loop  is  changed  to  a sine  wave,  the 
benefit  of  the  adaptive  Smith  Predictor  is  reduced,  as  shown  in  Figure  2-17.  Similarly  to 
Figure  2-11,  the  periodic  flowrate  Q(t)  profile  begins  at  time  2 hr,  has  a period  of  24  hr, 
and  has  an  amplitude  of  0.3  m /s. 

In  the  continuously  varying  flowrate  scenario,  only  the  linear  plant  simulation 
produces  a satisfactory  response.  As  in  Figure  2-11,  the  output  is  a smooth  sine  wave 
that  lags  behind  the  disturbance.  For  the  nonlinear  model,  the  controller  appears  to 
provide  even  poorer  control  than  with  a traditional  Smith  Predictor.  It  appears  poorer 
because  the  output  ranges  from  0.2  ppm  to  6.1  ppm.  However,  the  extremely  low  level  of 
0.2  ppm  is  only  observed  in  the  first  cycle.  In  the  second,  a minimum  of  0.5  ppm  is  again 
accomplished.  In  regard  to  reducing  the  size  of  the  peak  for  over-chlorination,  it  is 
slightly  lower  and  much  narrower,  both  of  which  are  improvements.  It  is  clear  that, 
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although  there  is  some  improvement,  a different  yet  controller  form  is  necessary  to 
achieve  satisfactory  performance. 

2.5  Novel  Process  Modification — Dynamic  Weir 

None  of  the  traditional  control  strategies  are  hilly  effective  when  the  flowrate  of 
influent  varies  throughout  the  range  typically  plaguing  wastewater  treatment  plants. 
Because  the  delay  is  so  large,  compensation  is  necessary,  yet  variations  in  the  delay  cause 
a Smith  Predictor  to  fail.  Primarily,  this  is  because  there  are,  in  fact,  two  variables  that 
must  be  controlled  to  guarantee  consistent  performance,  not  one.  The  first,  of  course,  is 
the  chlorine  concentration  itself,  and  the  second  is  the  residence  time  of  the  fluid  in  the 
reactor.  Because  the  fluid  is  experiencing  ongoing  chemical  reaction  for  as  long  as  it  is 
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retained  in  the  reactor,  it  is  necessary  to  control  how  long  it  remains  in  the  reactor  in 
order  to  guarantee  concentration  control. 

Controlling  the  residence  time  of  the  fluid  is  achieved  by  controlling  the  bulk 
velocity  of  the  fluid  while  in  the  reactor  (assuming  the  length  of  the  reactor  does  not 
change).  Recall  that  the  velocity  is  defined  as  v(t)  = Q(t)  / A (in  Section  2.2.3).  In  other 

words,  the  coefficient  of  — in  (2-28)  is  time-varying,  making  it  extremely  difficult  to 

oz 

predict  the  value  C(t,L)  for  random  v(t)  profiles.  However,  if  it  is  possible  to  manipulate 
A over  time  such  that  v(t)  = constant , significantly  better  control  of  chlorine 
concentration  may  be  achieved.  Physically,  adjusting  the  width  of  the  contact  basin  is  not 
likely  to  be  feasible;  the  height,  however,  could  be  easily  set  with  the  implementation  of  a 
movable  weir  located  at  the  end  of  the  reactor.  A reactor  with  such  a modification  is 
shown  in  Figure  2- 1 8 
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Figure  2-18.  Wastewater  chlorination  reactor  modified 
with  a dynamic  weir. 
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where  h(t)  is  the  height  of  the  fluid  in  the  reactor,  hx(t)  is  the  excess  height  of  fluid  above 

the  weir,  or  weir  head,  hw(t)  is  the  height  of  the  weir,  and  all  other  variables  are  the  same 
as  given  previously. 

It  is  assumed  that  the  height  of  fluid  in  the  reactor  is  the  same  along  the  entire 
length  of  the  reactor,  and  that  it  responds  uniformly  with  respect  to  z for  changes  in  either 
incoming  flowrate  or  weir  height.  The  assumption  is  based  upon  two  principles.  The 
first  is  the  Manning  equation  (Sturm,  2001),  which  calculates  the  slope  that  is  needed  to 
supply  the  driving  force  required  for  a given  flow 

c nV 

S = (2-57) 

where  S is  the  slope  of  the  channel  (dimensionless),  n is  a roughness  factor  specific  to  the 
channel  (dimensionless),  v is  the  velocity  of  the  water  in  the  channel  in  meters  per  second 
(m/s),  and  R is  the  hydraulic  radius  in  meters  (m).  Not  shown  is  a conversion  factor  of  1 
with  the  dimensions  necessary  to  make  the  equality  true.  Under  any  realistic  conditions, 
flow  of  such  low  velocity  will  only  cause  a height  gradient  of  a fraction  of  a centimeter 
per  meter  of  channel  length.  In  other  words,  at  steady-state,  h is  essentially  not  a function 
of  z. 

The  second  reason  that  h is  assumed  not  to  be  a function  of  z is  the  velocity  of 
propagation  of  small  disturbances  (Lencastre,  1987) 

u = Vgh  (2-58) 

where  u is  the  wave  velocity  in  meters  per  second  (m/s),  g is  the  gravitational  constant  in 

meters  per  second  squared  (m/  s2)  and  h is  the  depth  of  the  fluid  in  meters  (m).  Also 
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given  as  a form  of  the  Froude  number,  the  ratio  of  this  wave  velocity  u to  bulk  velocity  v 
shows  that  perturbations  in  the  system  pass  through  the  reactor  much  faster  than  the  fluid. 
In  other  words,  transients  attenuate  quickly,  and  the  system  can  be  assumed  to  be  of 
uniform  height  with  respect  to  location  in  the  reactor. 

2.5.1  Dynamic  Weir  Controller  Design 

2.5. 1.1  Embodiment  one — a model-based  ratio  controller 

As  implied  by  the  relationship  v(t)  = Q(t)  / A(t) , manipulation  of  height  is  simply 
a form  of  ratio  control.  Because  the  incoming  flowrate  is  already  being  measured,  the 
infonnation  is  readily  available.  It  is  next  necessary  to  consider  the  flowrate  of  fluid  over 
the  weir  at  the  exit  of  the  reactor.  There  are  many  different  flow  correlations  to  choose 
from,  and  each  contains  parameters  specific  to  each  plant.  In  practice,  significant 
investigation  may  be  required  to  determine  the  precise  flow  correlation  for  a given  plant. 
The  correlation  used  here  is  the  modified  Francis  weir  formula  (Geankoplis,  1993)  for  the 
case  of  no  suppressions  (i.e.  both  ends  of  weir  are  flush  with  the  side  of  the  reactor) 

Qw  =0.415whxl5V2g  (2-59) 

• 3 

where  Qw  is  the  flow  over  the  weir  in  cubic  meters  per  second  (m  /s),  g is  the 

2 

gravitational  constant  in  meters  per  second  squared  (m/s  ),  hx  is  the  weir  head  in  meters 
(m),  and  w is  the  width  of  the  basin  (m).  The  result  is  an  additional  state  for  the  system 

dV 

— = Q(t)-Qw  (2-60) 


which  upon  substitution  becomes 
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d7  = (Q(t) " °'4 1 5^w  ~ h*  (*»  h*  W 5 V2g) 


(2-61) 


In  addition,  we  must  recognize  that 


hx(t)  = h(t)  — hw  (t) 


(2-62) 


Substituting 


where  Q(t)  is  a measured  disturbance  and  hw  (t)  is  the  input.  By  modeling  the  height  as  a 
state,  it  is  now  possible  to  determine  if  transient  periods  where  Q Qw  affect  the 
behavior  of  the  reactor  significantly. 

The  final  obstacle  to  implementing  a ratio  control  is  knowledge  of  how  to  set  the 
height  of  the  weir  hw  (t)  based  upon  the  incoming  flowrate  Q(t)  to  yield  the  proper  fluid 

volume.  This  is  done  by  subtracting  a term  that  represents  the  height  of  the  head  for  that 
flowrate  from  the  total  height  desired. 


The  relationship  Q(hx)  is  given  in  (2-59);  however,  the  controller  requires  hx(Q) . 
Although  (2-59)  is  explicitly  invertible  (once  the  value  of  w is  known),  many  flow 
correlation  equations  are  not.  If  necessary,  the  inverse  may  be  approximated  with  a 
polynomial  fit;  for  example,  if  the  Kanapaha  facility  is  approximated  to  have  a width  of  8 


hw=h-hx(Q) 


(2-64) 


3 3 

m and  flowrate  in  the  range  0.1  m /s  < Q(t)  < 0.7  m /s,  the  weir  head  is  well-described  in 


terms  of  flowrate  by  the  equation 


hx(Q)  = — 0.0762Q(t)2  +0.2188Q(t)+.0156 


(2-65) 
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The  total  height  desired  must  also  be  calculated.  It  is  based  upon  the  desired  residence 
time,  x°es , which  set  by  the  operator  (i.e.  it  is  the  setpoint  of  the  ratio  controller).  It  is 
related  to  the  flowrate  and  the  volume  by 


T°  = 
res 


v(t) 

Q(t) 


where  V(t)  is  the  dynamic  volume.  Equivalently, 


tL  = Lw 


h(t) 

Q(t) 


Rearranging  and  substituting  into  (2-64)  yields 


(2-66) 


(2-67) 


hw(t)  = Q^C_h>(Q)  (2.68) 

where  hw  (t)  is  the  output  of  the  ratio  controller  in  meters  (m),  x°es  is  the  setpoint  of  the 
ratio  controller  in  seconds  (s),  and  hx(Q)  is  the  estimated  weir  head.  It  is  now  clear  that 
this  second  ratio  controller  is  in  fact  a residence  time  controller,  and  therefore  designated 


Figure  2-19.  Block  diagram  of  a wastewater  treatment 
plant  with  a dynamic  weir  and  model-based  ratio  controller 
under  cascade  chlorine  concentration  control. 
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RTC  (recall  that  the  input  concentration  is  also  implemented  through  ratio  control).  It  is 
incorporated  into  the  control  scheme  as  shown  in  Figure  2-19.  It  is  otherwise  identical  to 
the  previous  cascade  control  scheme.  The  RTC  also  uses  flowrate  Q(t)  as  an  input,  and 
calculates  the  desired  height  of  the  weir  hw  (t)  for  that  operating  condition.  This  signal  is 
the  passed  to  end  of  the  contact  basin,  where  it  is  implemented.  Finally,  as  in  the  case  of 
the  concentration  input  ratio  controller,  changes  in  the  weir  position  are  assumed  to  be 
fully  realized  far  faster  than  the  time  scale  of  the  system;  this  eliminates  the  need  to 
model  weir  movement  dynamics.  Note,  however,  that  the  total  height  of  the  fluid  is  now 
a state  of  the  model,  and  may  or  may  not  change  slowly  enough  to  affect  concentration. 

The  primary  advantage  of  this  embodiment  is  the  necessity  of  only  one 
measurement,  Q(t).  Coincidently,  that  measurement  is  frequently  already  available, 
producing  no  additional  burden  on  the  operator.  It  does,  however,  rely  upon  a model  of 
the  flow  over  a weir.  Since  the  physical  dimensions  of  each  wastewater  treatment  plant 
are  unique,  the  flow  correlation  must  be  identified  on  a case-by-case  basis.  Exhaustive 
experimentation  may  be  required  to  define  it.  In  addition,  if  the  model  is  not  accurate, 
degradation  of  controller  performance  may  be  may  be  observed. 

2.5.1.2  Embodiment  two — a feedback  PI  controller 

A standard  PI  controller  can  also  be  utilized  to  achieve  residence  time  control.  In 
digital  position  form,  it  is  given  as 

K (0  = h|T  + Kc[e(t)  + — £e(t  - kT)]  (2-69) 

TI  k=0 

where  Kc  is  the  proportional  constant  of  the  controller,  x,  is  the  integral  time  constant  of 

the  controller,  h^is  the  reference  position  of  the  weir,  and  T is  the  sampling  period. 
Furthermore,  the  feedback  error  e(t)  is  calculated  as  follows 
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e(t)  = C(t)-LW-^-  (2-70) 

where  e(t)  is  interpreted  as  the  difference  between  the  set-point  x°es  and  the  estimated 
residence  time  xres  given  as  the  last  term  in  (2-70).  Again  it  clear  that  controller  is 


Figure  2-20.  Block  diagram  of  a wastewater  treatment 
plant  with  a dynamic  weir  and  Feedback  PI  controller 
under  cascade  chlorine  concentration  control. 

charged  with  residence  time  control,  and  this  embodiment  is  also  referred  to  as  a RTC.  It 
is  incorporated  into  the  control  scheme  as  shown  in  Figure  2-20,  which  is  identical  to 
Figure  2-19  except  for  the  addition  of  a fluid  height  measurement. 

In  this  formulation,  the  RTC  again  uses  flowrate  Q(t)  as  an  input,  and  additionally 
uses  fluid  height  h(t)  as  a second  input.  With  these  two  pieces  of  information,  it 
calculates  the  current  residence  time  for  fluid  in  the  reactor.  The  PI  controller  tries  to 
then  eliminate  any  discrepancy  between  this  value  and  the  setpoint  by  sending  a signal  to 
adjust  the  height  of  the  weir.  As  in  the  case  of  the  model-based  ratio  controller,  changes 
in  the  weir  position  are  assumed  to  be  fully  realized  far  faster  than  the  time  scale  of  the 
system;  this  eliminates  the  need  to  model  weir  movement  dynamics.  Finally,  the  total 
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height  of  the  fluid  is  again  now  a state  of  the  model,  and  may  or  may  not  change  slowly 
enough  to  affect  concentration. 

The  primary  advantage  of  this  embodiment  is  that  no  flow  correlation  is  needed; 
in  contrast,  recall  that  the  model-based  ratio  controller  (i.e.  Embodiment  1)  requires  one. 
However,  in  addition  to  a measurement  of  Q(t),  it  needs  a measured  or  estimated  value  of 
h(t) — something  not  needed  by  the  ratio  controller.  Furthermore,  the  additional  PI 
controller  will  require  tuning.  Again,  since  the  physical  dimensions  of  each  wastewater 
treatment  plant  are  unique,  exhaustive  experimentation  may  be  required  to  achieve  the 
proper  parameters. 

Finally,  for  either  embodiment,  it  must  be  decided  what  the  bounds  of  weir  height 
are.  The  upper  bound  is  modeled  simply  as  the  present  height  of  the  fixed  weir.  This  is 
because  it  can  be  assumed  that  plants  are  sized  according  to  their  maximum  expected 
capacity.  For  that  maximum  value,  they  must  be  large  enough  to  retain  the  fluid  some 
minimum  amount  of  time.  Therefore,  it  is  never  expected  that  the  flow  will  go  higher, 
causing  the  prescribed  weir  height  to  exceed  the  present  value.  A very  different 
consideration  governs  the  choice  for  minimum  weir  height.  It  is  limited  either  by 
physical  on-site  constraints,  or  by  a limit  of  how  shallow  a channel  of  water  may  become 
before  uniform  flow  assumptions  break  down.  It  has  been  suggested  in  literature  that  a 
realistic  lower  bound  is  1 ft,  so  0.3  m is  chosen  for  the  purposes  of  simulation. 

2.5.2  Closed-Loop  Response  of  the  Entire  Model  for  Step  Changes  in  Flowrate 

Using  a Dynamic  Weir 

By  reusing  the  tuning  parameters  given  for  a cascade  chlorine  concentration 
control  system  with  standard  PI  controllers  (Table  2-8),  any  improvement  gained  by 
implementing  a dynamic  weir  can  be  observed.  In  the  case  of  the  model-based  ratio 
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controller,  this  constitutes  all  necessary  controller  parameters.  However,  for  the  PI 
feedback  system,  the  RTC  itself  is  an  additional  PI  controller,  and  requires  tuning;  trial- 
and-error  tuning  methods  show  that  Kw  =0.01  m/s  and  xw  = 1000s  are  effective 
choices  for  these  parameters  for  this  plant.  Figure  2-21  records  the  response  of  the 
system  equipped  with  a model-based  ratio  weir  controller  when  subjected  to  the  same 
series  of  step  changes  described  in  the  previous  section:  a small  early  step-increase  in 
disturbance  Q(t)  at  time  2 hr,  then  a large  step-decrease  in  disturbance  at  time  12  hr,  and 
finally  another  step-increase  at  time  34  hr,  returning  to  the  original  flowrate  Q(t).  The 
first  plot  shows  the  output  concentration  (solid  line)  and  the  setpoint  (dotted  line)  as 
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Figure  2-21.  Closed-loop  response  of  model  with  a pair  of 
chlorine  concentration  controllers  arranged  in  cascade 
configuration  and  a model-based  ratio-controlled  dynamic 
weir,  subjected  to  a step  changes  in  flowrate. 
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functions  of  time.  A value  of  2 ppm  is  again  chosen  for  the  setpoint,  with  1 ppm  being 
the  minimum  allowable  threshold.  The  second  plot  shows  the  chlorine  concentration 

entering  the  contact  area  Co(t)  for  each  of  the  two  models  over  time.  The  third  plot  is 

new,  and  gives  the  height  of  the  weir  as  a function  of  time.  Finally,  the  disturbance  Q(t) 
is  shown  in  the  fourth  plot  as  a function  of  time.  In  addition,  the  setpoint  of  the  RTC  is 
chosen  as  x°es  = 3600  s . 

Figure  2-22  is  the  identical  simulation  performed  for  the  case  when  the  RTC  is  a 
PI  feedback  controller.  In  other  words,  a value  of  2 ppm  is  chosen  for  the  setpoint,  with  1 
ppm  being  the  minimum  allowed  value.  Also,  the  setpoint  of  the  RTC  is  x°es  = 3600  s . 

Finally,  each  of  the  plots  in  Figure  2-22  represents  the  same  value  it  did  as  explained  for 
Figure  2-21. 

The  responses  of  the  two  types  of  RTCs  are  almost  identical.  Furthermore,  it  is 
clear  that  the  dynamic  weir  arrangement  yields  an  excellent  chlorine  concentration 
control  response  for  the  conditions  described.  They  do  an  excellent  job  returning  the 
system  to  setpoint,  only  deviating  approximately  0.2  ppm  in  the  most  extreme  case. 
Furthermore,  they  achieve  the  disturbance  rejection  at  most  in  about  three  hours,  much 
faster  than  the  case  with  a fixed  weir  height.  This  behavior  confirms  that  the  dynamics  of 
the  system  are  fast  enough  to  ensure  that  gravity  flow  should  be  sufficient  to  deliver  the 
system  to  the  desired  height  in  a reasonable  amount  of  time. 

Interestingly,  the  most  problematic  step  change  is  now  a large  step-up,  in  contrast 
to  the  fixed-weir  case,  where  a dramatic  step-down  causes  severe  control  problems.  This 
cause  of  this  phenomenon  is  the  uncontrolled  nature  of  the  inlet  flowrate.  When  flowrate 
surges  from  a low-flow  to  a high-flow  condition,  and  the  weir  attempts  to  quickly  create  a 
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Figure  2-22.  Closed-loop  response  of  model  with  a pair  of 
chlorine  concentration  controllers  arranged  in  cascade 
configuration  and  a PI  feedback-controlled  dynamic  weir, 
subjected  to  a step  changes  in  flowrate. 

greater  volume  of  fluid  in  the  reactor,  the  rate  at  which  it  can  be  accomplished  is  fully 
dependent  upon  how  quickly  the  new  flowrate  can  supply  the  needed  amount  of  fluid. 
This  fill-time  can  be  approximated  using  the  formula 


tfiii  = 


Lw(hf-hj) 

Qr 


(2-71) 


where  hf  and  hf  are  the  final  and  initial  heights,  respectively,  and  Qf  is  the  final 
flowrate.  The  maximum  fill-time  is  specific  each  plant,  depending  upon  the  values  of  the 
parameters  given  in  (2-71).  As  an  example,  for  the  plant  detailed  throughout  this  study, 
which  is  described  by  L = 1 20  m , w = 8 m , hf  = 2.62  m , h;  = 0.37  m , and 
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Qf  - 0.7  m3  / s,  it  is  calculated  that  tfl]1  = 51  minutes.  While  this  is  a highly  undesirable 
situation,  it  should  be  noted  that  such  a surge  in  flowrate  is  equally  undesirable  in  a plant 
with  a fixed  weir,  and  even  when  added  to  the  described  residence  time  (60  minutes),  the 
resulting  fluid  is  retained  for  less  time  than  is  frequently  observed  in  plants  with  a fixed 
weir. 

2.5.3  Closed-Loop  Response  of  the  Entire  Model  for  Varying  Flowrate  Using  a 
Dynamic  Weir 

Figure  2-23  shows  the  performance  of  the  standard  cascade  control  system  (the 
tuning  parameters  remain  unchanged)  supplemented  with  a feedforward  RTC  and  a 


t [hrs] 


Figure  2-23.  Closed-loop  response  of  the  plant  with  a pair 
of  chlorine  concentration  controllers  arranged  in  cascade 
configuration  and  a model-based  ratio-controlled  dynamic 
weir,  subjected  to  varying  flow. 
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dynamic  weir  when  the  flowrate  is  modeled  as  the  previously  described  sinusoidal  wave. 
As  in  previous  sections,  the  sine  wave  begins  at  time  2 hr,  has  a period  of  24  hr,  and  has 

an  amplitude  that  is  three-fourths  of  the  nominal  flowrate  (0.4  m3/s),  or  0.3  m3/s.  The 

3 3 

result  is  again  a disturbance  Q(t)  that  varies  from  0.1  m /s  to  0.7  m /s.  As  in  the  previous 

section,  the  setpoint  of  the  second  ratio  controller  is  chosen  as  x°es  = 3600  s . 

Similarly  to  Section  2.5.2,  Figure  2-24  is  exactly  the  same  simulation  performed 
to  generate  Figure  2-23,  except  that  the  RTC  is  instead  a PI  feedback  controller  with 
tuning  parameters  Kw  =0.01  m/s  and  rw  = 1000s. 


t [hrs] 


Figure  2-24.  Closed-loop  response  of  the  plant  with  a pair 
of  chlorine  concentration  controllers  arranged  in  cascade 
configuration  and  a PI  feedback-controlled  dynamic  weir, 
subjected  to  varying  flow. 
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Again,  it  is  clear  that  the  performance  of  the  two  different  RTC  embodiments  is 
practically  identical.  In  comparing  this  result  to  other  control  strategies  described,  it  is 
clear  that  the  value  of  CL  is  controlled  much  more  tightly,  falling  in  the  approximate 
range  1 .99  ppm  < CL  (t)  < 2.0 1 ppm . This  is  in  contrast  0.5  ppm  < (t)  < 6.0  ppm  in 

the  simulation  of  a weir  with  fixed  height  shown  in  Figure  2-13.  This  means  that  a 
dynamic  weir  can  virtually  eliminate  the  effect  that  flowrate  fluctuations  have  on  the 
outlet  chlorine  concentration. 

It  is  important  to  note  that  perfect  control  cannot  be  expected.  Only  the  most 
dramatic  disturbance  (variable  flow  and  its  effect  upon  residence  time)  has  been  modeled; 
in  practice,  the  system  suffers  from  other,  less  severe  disturbances.  The  simulation  study 
serves  to  show  that,  even  while  ignoring  lesser  disturbances,  a traditional  control  strategy 
is  inadequate  to  control  a system  with  a fixed  weir.  The  dynamic  weir  solves  problem  of 
the  most  severe  disturbance  (variable  flow,  and  therefore,  residence  time),  leaving  only 
less  severe  disturbance-related  problems  to  be  addressed. 

2.5.4  Robustness  Analysis  for  the  Model-Based  Residence  Time  Controller 

It  should  be  noted  that  no  weir  correlation  is  exactly  accurate.  Similar  to 
measurement  error,  when  the  actual  flow  over  the  weir  differs  from  that  which  is 
predicted  by  the  correlation,  performance  may  be  degraded.  In  order  to  quantify  how  far 
the  calculated  control  variable  (i.e.  the  residence  time)  can  be  from  its  true  value, 
consider  error  possible  in  the  model  of  flow  over  a weir  as  a function  of  weir  head  h x (t) . 
Rearranging  (2-62)  and  substituting  into  (2-67) 


Lw 


hw(t)  + hx(t) 


Q(t) 


(2-72) 
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If  the  weir  correlation  being  used  (2-59)  can  be  inverted  directly,  the  error  can  be 
cast  as  a multiplicative  factor  of  the  correlation 


Qw  = 0.415w  hx‘  5^/2g(l  + s)  (2-73) 

where  Qw  is  the  actual,  or  true,  flowrate  and  lOOe  is  the  percent  error  in  weir-flow  rate. 
Rearranging  this  equation  to  yield  hx  in  terms  of  Qw , the  relationship  is 


h'x  = (: 


Qv 


0.415w^/2g 

Equivalently,  hx  can  be  represented  in  terms  of  h and  e 


2 2 
)3(l  + e)  3 


(2-74) 


_2 

hx=hx(l  + £)3  (2-75) 

Substituting  (2-75)  into  (2-72)  to  find  the  true  residence  time  yields  the  equation  used  for 

control  design  at  steady-state 


Tres(£)  = 


Lw 

QCt) 


2 \ 


hw(t)  + hx(t)(l  + e)  3 


(2-76) 


Next,  solve  (2-72)  for  hw  (t)  in  the  case  where  xres(e  = 0)  (i.e.,  xres  = x°es) 


hw(t)  = -'7Q(t)-hx(t)  (2-77) 

Lw 

and  substituting  into  (2-76),  after  rearrangement,  yields  xres(e)  in  terms  of  additive 
perturbations  in  x^ 


X (e)  = x°  +i^M> 
' "*  Q(t) 


( 


2 ^ 


(1  + e)  3 -1 


(2-78) 


which  can  also  be  expressed  in  terms  of  percent  error  from  the  target  value 
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ATres(e)  = 


Lwhx(t) 

x?«Q(0 


_2  \ 

(l  + £)  3 -1 


(2-79) 


where 


Atres(£):=  (xres(s)  - tL)  / . (2-80) 

Finally,  substituting  (2-67)  into  (2-79)  and  canceling  terms  yields  the  expression 


Axres(£)  = 


MO 


h(t) 


2 \ 


(1  + e)  3-l 


(2-81) 


Although  typical  weir  correlations  are  accurate  to  within  a few  percent,  the 
robustness  of  this  control  design  can  be  shown  by  plotting  wide  range  of  errors.  Figure 


Figure  2-25.  Plot  of  observed  per  cent  error  in  xres  for 
various  values  of  weir  correlation  error  e . 
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2 25  illustrates  the  results  for  low,  middle,  and  high-flow  conditions.  In  other  words,  the 
flowrate  Q(t)  is  fixed  at  three  separate  points,  and  a curve  is  generated  for  each  by 
varying  the  value  of  £ . Each  calculation  assumes  that  L = 120  m,  w = 8 m,  and 

^res  3600s . In  addition,  the  weir  height  hw(t)  is  set  at  the  appropriate  height  (but  not 
reported)  for  each  of  the  three  flowrates. 

It  is  clear  that  the  controller  is  extremely  robust,  given  that  it  can  easily  tolerate 
errors  in  the  weir  correlation  of  up  to  50%  while  yields  errors  in  residence  time  of  < 6%. 
The  reason  is  that  in  general,  h(t)  » hx(t) , so  errors  in  hx(t)  simply  do  not  have  an 
extremely  significant  effect.  The  trend  is  evident  in  (2-81).  This,  however,  underscores 
the  importance  of  having  an  accurate  value  of  hw  (t) ; therefore,  it  is  critical  to  monitor 
whether  effects  such  as  hysteresis  cause  significant  error  in  weir  position.  If  the  are  any 
problems,  however,  the  problem  is  easily  remedied  by  making  the  RTC  control  signal  the 
setpoint  to  a slave  PI  controller  designed  to  guarantee  that  the  weir  is  in  exactly  the 
position  prescribed  by  the  RTC.  Finally,  while  the  robustness  analysis  would  be  more 
complicated  it  the  correlation  is  not  invertible,  it  can  be  assumed  from  the  insensitivity  of 
the  system  to  correlation  errors  that  the  result  would  not  be  significantly  different. 

2.5.5  Propagation  of  Error  Analysis  for  the  PI  Feedback-Controlled  System 

Errors  in  the  measurement  ol  key  variables  can  cause  degradation  of  controller 
performance.  In  the  model-based  embodiment  of  the  controller,  the  analysis  is  relatively 
simply.  The  true  value  of  the  only  measurement,  Q(t),  can  be  characterized  to  fall  within 
certain  bounds  above  or  below  the  observed  value.  These  high  and  low  values  directly 
limit  the  range  of  possible  true  residence  time  for  a given  moment  in  time,  since  all  other 
variables  are  assumed  to  be  known  exactly. 
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In  contrast,  the  PI  feedback  controller  requires  two  measurements,  and  each 
contributes  to  the  overall  error  separately.  In  addition,  if  the  two  variables  are  not 
independent,  then  a convoluted  interaction  can  contribute  as  well.  The  possible  range  of 
error  is  therefore  significantly  more  complicated  than  in  the  single  variable  case.  One 
way  of  quantifying  it  is  through  a propagation  of  variance  analysis,  which  is  given  by  the 
equation 


= 


SG"| 

dQj 


fd  gn2 


°*Q+I 


Q,h 


V dh  y 


CTh  + 


Q.h 


( 

V^Qy 


Q.h 


fdG 
dh  J 


CT 


Oh 


(2-82) 


Q.h 


is  a 


where  ctt  is  the  variance  of  residence  time,  aQ  is  the  variance  of  the  flowrate 
measurement,  ah  is  the  variance  of  the  height  measurement,  aQh  is  the  covariance 

between  the  measurements  of  Q and  h,  Q and  h are  operating  point  values,  and  G is 
function  that  describes  the  relationship  between  the  variables.  Determining  the  values  of 
<Jq  and  ah  are  straightforward,  and  achieved  by  considering  the  specifications  given  for 
the  appropriate  measurement  device. 

For  a representative  magnetic  flow  meter  of  appropriate  capacity,  the 
specification  is  given  as  ( ± 0.5  % of  flowrate).  Consequently,  the  error  observed  for  a 
measurement  with  this  device  is  a function  of  operating  point  ( i.e . the  value  of  Q). 
Assuming  that  the  observed  values  from  the  flow  meter  fall  within  three  standard 
deviations  of  the  true  value 


3ctq  = 0.005Q 

which  can  easily  be  substituted  into  (2-82)  for  different  values  of  Q . 


(2-83) 
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In  examining  the  variance  of  the  height  measurement  crh , it  is  found  that  error  of 
a representative  piece  of  hardware  for  the  appropriate  height  range  is  given  as 
(±0.1  % of  span).  In  contrast  to  the  flow  criterion,  this  value  does  not  change  for 
different  operating  conditions,  since  it  is  assumed  that  the  height  measurement  device 
will  be  subject  to  constant  maximum  and  minimum  values,  regardless  of  the  present 
value  of  the  variable.  Therefore,  in  the  case  of  the  described  plant,  with  the 
aforementioned  maximum  and  minimum  height  values  of  0.3  and  3.0  m,  respectively, 
and  the  same  standard  deviation  assumption  given  for  the  flow  meter,  it  found  that 

3a^=  0.0027m  (2-84) 

which  can  also  be  subsequently  substituted  into  (2-82). 

The  identification  of  the  value  of  ctQIi  is  generally  more  difficult  that  either  oQ 
or  ah . It  is  necessary  to  investigate  how  one  affects  the  other.  In  the  case  of  the 
wastewater  treatment  plant,  the  value  of  Q directly  affects  the  value  of  h,  and  their 
covariance  agh  can  be  estimated  by  generating  Q and  h signals  of  appropriate  variance 

(as  defined  by  (2-83)  and  (2-84))  and  sufficient  length  (vectors  of  1000  points  were 

Table  2-11.  Summary  of  the 
square  of  the  covariance  for  low, 
middle,  and  high-flow  operating 
points. 


Q (m3/s) 

2 

aQh 

0.1 

4.719E-9 

0.4 

2.132E-8 

0.7 

3.927E-8 
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generated  in  Matlab  for  this  exercise)  and  calculating  the  resulting  covariance.  When 
performed  for  the  previously  described  low,  middle,  and  high-flow  conditions  (0. 1 , 0.4, 
and  0.7  m3/s),  the  observed  standard  deviations  are  as  given  in  Table  2-1 1. 

Finally,  the  function  G must  be  defined.  It  is  simply  the  relationship  between  the 
variables,  which  is  in  reality  a slightly  more  general  form  than  is  already  given  by  (2-67) 


x 


res 


= Lw 


h(t) 

Q(t) 


(2-85) 


Partial  derivatives  of  (2-85)  produce  the  remaining  terms  in  (2-82),  with  the  appropriate 


Table  2- 1 2.  Summary  of  the  range 
of  error  in  xres  for  low,  middle, 
and  high-flow  operating  points. 


Q (m3/s) 

h (m) 

Error  in  xres  j 

0.1 

0.337 

± 31  s 

0.4 

1.409 

± 19s 

0.7 

2.492 

± 18  s 

values  of  Q and  h used  for  each  operating  condition,  as  given  in  Table  2-12.  Finally, 
the  previously  given  values  of  L =120m  and  w = 8 m are  also  used.  Considering  that 
the  setpoint  for  these  simulations  is  3600  s,  it  is  clear  that  the  measurement  error  tends 
introduce  a deviation  of  less  than  a percent. 

2.5.6  Additional  Benefits  of  the  Dynamic  Weir  Design 

There  are  benefits  in  addition  to  the  highly  improved  performance  associated  with 
using  a dynamic  weir  to  control  residual  chlorine  concentration.  Primarily,  it  is  much 
easier  to  guarantee  sufficient  levels  of  microbe  disactivation,  thereby  virtually 
eliminating  the  risk  of  exposing  humans  to  numerous  varieties  of  bacteria,  viruses  and 
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protozoa,  including  Salmonella,  Shigella  and  Vibrio  cholera.  Furthermore,  with  chlorine 
levels  no  longer  reaching  the  5 to  10  ppm  level,  or  higher,  the  risk  of  THM  formation  is 
greatly  reduced  as  well. 

Equally  attractive  to  the  wastewater  treatment  industry  is  the  potential  of  the 
dynamic  weir  to  significantly  reduce  operating  costs.  Some  of  the  reduction  can  be 
realized  indirectly;  for  instance,  the  need  for  equalization  basins  (Foess  et  al.,  1976; 
Ongerth,  1979),  which  are  designed  to  attenuate  the  highly  varied  flow  into  the 
chlorination  basin,  is  gone.  In  addition,  it  should  be  much  easier  to  avoid  fines  for  non- 
compliance.  Furthermore,  if  sophisticated,  expensive  advanced  control  techniques  are 
used  to  control  a particular  plant,  they  can  be  abandoned  in  favor  of  the  dynamic  weir 
with  its  easily  programmed  residence  time  controller. 

However,  the  most  notable  economic  benefit  is  that  less  chlorine  is  required  to 
achieve  the  same  level  of  safety.  (Of  course,  the  dynamic  weir  can  be  used  with  any 
chemical  treatment;  chlorine  is  simply  used  throughout  this  study  for  both  simplicity  and 
consistency.)  For  example,  assuming  a constant  ammonia  chlorine  demand  of  3 ppm,  the 
four  basin  designs  covered  previously  require  a daily  supply  of  chlorine  as  given  in  Table 
2-13.  Each  simulation  was  again  run  for  50  hours,  with  the  sinusoidal  flow  beginning  at  t 
= 2 hrs,  and  a day-long  period  with  flow  ranging  between  0.1  and  0.7  m3/s.  This 
corresponds  approximately  to  an  8 mega-gallon  per  day  (MGD)  average  plant  that  is 
subjected  to  fluctuations  from  to  2 to  14  MGD,  respectively.  Total  chlorine  use  for  a 
single  day  is  calculated  simply  by  integrating  the  dosing  flowrate  for  the  second  24  hr 
flow  cycle;  as  before,  this  is  done  to  eliminate  any  transient  effects  caused  by  beginning 
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the  simulation  at  steady-state  with  respect  to  flowrate.  Finally,  it  is  assumed  that  the  cost 
of  chlorine  gas  is  $320.00  per  US  ton. 


Table  2-13.  Comparison  of  the  daily  chlorine  usage  cost 
for  each  of  the  plant  configurations  studied. 


Weir  Type 

fixed 

fixed 

fixed 

dynamic 

Master  Controller 
Type 

PI 

PI  w/  Smith 
Predictor 

PI  w/  Adaptive 
SP 

PI 

1 Daily  Cost  ($) 

100 

98 

93 
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The  observed  chlorine  usage  reduction  is  therefore  approximately  25%.  The 


percentage  savings  can  vary,  though,  if  the  actual  ammonia  chlorine  demand  deviates 
from  3 ppm.  However,  the  absolute  savings  remain  the  same.  For  a plant  of  the  given 
size,  the  dynamic  weir  system  will  save  the  treatment  plant  approximately  $10,000  per 
year  over  the  most  common  type  of  system  utilizing  a fixed  weir  (i.e.  one  that  uses  a PI 
master  controller).  Furthermore,  with  much  more  effective  control  possible,  the  setpoint 
{i.e.  presently  2 ppm)  may  be  moved  closer  to  the  allowable  threshold  (i.e.  1 ppm  in  this 
case),  thereby  saving  even  more  in  chlorine  costs  (perhaps  an  additional  $5000  yearly) 
without  the  risk  of  failure  to  comply  with  regulations.  This  savings  is  realizable 
primarily  because,  even  when  flow  is  equalized  in  a fixed-volume  reactor,  the  reactor  is 
still  not  being  operated  optimally  in  the  sense  that  the  residence  time  of  wastewater 
within  the  reactor  is  greater  than  it  needs  to  be,  requiring  more  chlorine  for  disinfection. 

Finally,  reduced  chlorine  use  alone  provides  safety  benefits.  With  the  demand  for 
chlorine  reduced,  it  is  necessary  to  store  less  on-site,  which  in  turn  would  likely  reduce 
the  scale  of  the  problem  should  an  accident  occur. 
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2.6  Conclusions 

The  effect  of  ammonia  nitrogen  demand  dictates  that  failure  to  utilize  cascade 
control  will  most  likely  result  in  degraded  outlet  chlorine  concentration  performance. 
Hence,  a cascade  controller  is  needed;  however,  for  dosage  adjustment,  a cascade  control 
strategy  alone  does  not  guarantee  acceptable  performance,  primarily  because  it  does  not 
address  the  effect  that  varying  flowrate  has  on  residence  time.  The  implementation  of  a 
residence  time  controller  (RTC)  supplementary  to  existing  dosage  control  techniques  can 
bring  about  excellent  perfonnance. 

Because  the  RTC  scheme  reduces  the  observed  variation  in  residence  time  from 
several-fold  to  just  a few  percent,  the  pre-existing  control  structure  becomes  much  more 
effective  in  delivering  acceptable  outlet  chlorine  concentration  control.  A fundamental 
explanation  for  this  trend  is  that  a disturbance  in  the  form  of  varying  flowrate  has  two 
distinct  and  independent  effects  on  the  system.  The  first  is  the  interruption  of  constant 
input  concentration,  an  undesirable  effect  that  is  commonly  rejected  using  a ratio 
controller.  However,  the  flowrate  fluctuations  also  cause  a second  effect — varying 
residence  time,  for  which  the  concentration  ratio  controller  is  ineffective.  The 
feedforward  RTC  can  be  interpreted  as  a novel  form  of  a ratio  control,  with  the  task  of 
maintaining  approximately  constant  the  ratio  of  fluid  height  in  the  reactor  to  flowrate  (2- 
67),  resulting  in  an  approximately  constant  residence  time  of  fluid  in  the  reactor.  Thus, 
although  there  are  two  primary  disturbances,  namely  Q(t)  and  5(t),  they  induce  three 
effects  on  the  system  ( i.e .,  (i)  input  concentration  interruption,  (ii)  variable  residence 
time,  and  (iii)  a load  on  the  system  that  consumes  chlorine  in  the  dosage  area);  therefore, 
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three  controllers  (ratio,  RTC,  and  slave  controllers)  are  required  to  adequately  reject  all 
effects  of  the  disturbances. 

Both  the  feedforward  and  feedback  configurations  of  the  RTC  perform  adequately 
when  the  measured  data  is  not  filtered.  It  is  anticipated,  however,  that  performance  can 
be  improved  by  employing  a digital  filter.  In  fact,  a first-order  exponential  filter 
(Ogunnaike  and  Ray,  1994)  with  filter  parameter  A = 0.01  should  yield  satisfactory 
performance.  Despite  the  low  value  of  the  parameter,  data  can  be  taken  several  times  per 
second;  high-frequency  data  acquisition  assures  that  no  significant  additional  delay  is 
introduced  by  the  filter,  since  the  residence  time  requirement  is  likely  to  be  thirty  minutes 
or  more. 

There  are  two  alternative  ways  to  implement  the  RTC  scheme:  namely,  a 
feedforward  model-based;  and  a feedback  non-model-based  configuration.  Trade-off 
exists  between  them,  and  plant-specific  considerations  should  be  used  to  determine  the 
best  control  configuration.  For  example,  the  feedforward  algorithm  requires  an  inverted 
model  of  flow  over  a weir.  In  contrast,  the  feedback  algorithm  does  not  require  a model, 
but  does  require  a fluid  height  measurement  as  well  as  controller  tuning.  The  ability  of 
the  engineer  to  acquire  the  requisite  information  ultimately  determines  the  appropriate 
choice. 

Both  the  feedforward  and  feedback  RTC  schemes  are  capable  of  yielding 
excellent  control,  provided  that  the  weir  is  placed  in  the  proper  position.  It  is  possible 
that  the  weir  will  not  always  move  to  the  position  prescribed  by  an  RTC;  however,  the 
problem  is  easily  remedied  by  making  the  RTC  control  signal  become  the  setpoint  to  a 
slave  PI  controller  designed  to  guarantee  that  the  weir  is  in  exactly  the  position  prescribed 
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by  the  RTC.  Furthermore,  the  feedforward  RTC  is  an  extremely  robust  controller  with 
respect  to  weir-flow  modeling  error,  making  weir  head  estimation  much  less  critical  than 
weir  position  in  terms  of  guaranteeing  good  performance.  Finally,  while  the  robustness 
analysis  would  be  more  complicated  if  the  correlation  for  flow  over  a weir  is  not 
explicitly  invertible,  it  can  be  expected  from  the  insensitivity  of  the  system  to  correlation 
errors  that  the  result  would  not  be  significantly  different. 

The  most  notable  economic  benefit  is  that,  under  the  conditions  of  the  simulation, 
the  RTC  scheme  has  the  potential  to  realize  an  average  25-30%  annual  reduction  in  the 
amount  of  chlorine  used  while  facilitating  better  control.  This  is  possible  because 
presently,  virtually  all  wastewater  passing  through  the  disinfection  tank  remains  far 
longer  than  required  by  the  regulations;  consequently,  the  water  must  be  dosed  with  more 
chlorine  than  would  otherwise  be  necessary.  For  example,  assuming  an  average 
ammonia  nitrogen  chlorine  demand  5(t)  of  3 ppm,  the  fixed-weir  cascade  design  with  a 
PI  master  controller  and  moving  weir  with  a feedforward  RTC  control  design  are 
compared,  and  found  to  require  an  annual  supply  of  chlorine  of  1 15  and  83  US  tons  , 
respectively  (i.e.,  a difference  of  32  US  tons).  To  determine  yearly  usage,  simulation  of 
each  control  shceme  is  again  run  for  50  hours,  with  the  sinusoidal  flow  beginning  at  t = 2 
hrs,  and  a day-long  period  with  flow  ranging  between  0.1  and  0.7  m3/s.  This  corresponds 
approximately  to  an  8 mega-gallon  per  day  (MGD)  average  plant  that  is  subjected  to 
fluctuations  from  to  2 to  14  MGD.  Total  daily  chlorine  use  is  calculated  by  integrating 
the  dosing  flowrate  for  the  second  24  hr  flow  cycle;  this  is  done  to  eliminate  any  initial 
transient  effects.  The  daily  consumption  is  then  extrapolated  to  a yearly  usage. 
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The  extent  of  percentage  chlorine  savings  can  vary,  however,  if  the  average 
ammonia  chlorine  demand  deviates  from  3 ppm.  Nonetheless,  the  absolute  savings 
remain  the  same.  In  other  words,  the  amount  of  chlorine  presently  required  to  react  with 
the  ammonia  nitrogen  in  the  dosage  area  must  remain  the  same  under  all  control  schemes; 
however,  the  required  amount  of  chlorine  remaining  available  in  the  system  for 
disinfection  reactions  in  the  contact  area  is  significantly  reduced.  For  a plant  of  the  given 
size,  the  dynamic  weir  system  can  save  the  treatment  plant  approximately  32  US  tons  per 
year  over  the  fixed-weir  design.  Furthermore,  with  much  more  effective  control  possible, 
the  setpoint  (/.<?.,  presently  2 ppm)  may  be  moved  closer  to  the  allowable  threshold  (i.e.,  1 
ppm  in  this  case),  thereby  saving  even  more  in  chlorine  costs  without  the  risk  of  failure  to 
comply  with  regulations.  Finally,  while  the  setpoint  of  the  RTC  is  3600  s throughout  this 
study,  regulations  are  usually  closer  to  1800  s,  providing  for  the  potential  to  save  even 
more  chlorine.  For  example,  additional  simulations  show  that  the  yearly  average  chlorine 
usage  for  a plant  this  size  with  a new  outlet  concentration  setpoint  of  1.2  ppm  and  a 
residence  time  setpoint  of  1800  s would  be  further  reduced  to  approximately  63  US  tons. 
The  yearly  total  chlorine  savings  for  an  8 MGD  plant  then  becomes  approximately  52  US 
tons,  which  represents  a savings  ceiling  of  up  to  45%. 


CHAPTER  3 

STABILITY  OF  PI-CONTROLLED  RUN-TO-RUN 
STATIC  SYSTEMS  WITH  SIGNIFICANT  DELAY 

3.1  Introduction 

There  are  numerous  industrial  processes  that  may  be  automated  using  run-to-run 
control  schemes.  Such  processes  are  inherently  discrete,  and  consequently,  they  are  very 
similar  to  batch  systems,  with  the  subtle,  yet  important,  distinction  that  information  about 
each  batch  is  only  measurable  after  the  batch  is  completed.  The  semiconductor  industry 
is  rich  with  examples  of  run-to-run  control  applications,  such  as  plasma  etching  (Card  el 
al.,  1998)  and  polishing  (Boning  et  al.,  1996;  Moyne  et  al.,  2001).  Additional  examples 
include  the  serial  pressing,  sputtering,  and  coating  operations  in  compact  disc  production 
(Gillet  el  al.,  2001)  and  the  roasting,  cooling,  and  screening  stages  that  are  necessary  for 
peanut  and  hazelnut  processing  (Nagaraju  el  al.,  1997;  Ozdemir  and  Devres,  2000; 
Ozdemir  and  Devres,  2000). 

In  addition,  many  unit  operations  that  are  automated  with  run-to-run  control  are 
contained  within  of  a series  of  production  stages.  Therefore,  a significant  characteristic 
of  the  run-to-run  control  schemes  considered  here  is  that  properties  created  in  an  early 
stage  cannot  be  measured  immediately  after  that  stage;  rather,  they  must  complete  a 
series  of  steps  before  data  can  be  gathered  (Kosut  et  al.,  1998).  When  coupled  with  large 
throughputs,  a large  measurement  delay  may  result.  Additionally,  the  assumption  is 
made  in  this  study  that  dynamic  effects  from  one  batch  to  the  next  are  negligible.  In 
reality,  some  limited  dynamic  effects  may  linger  from  one  run  to  the  next,  such  as  latent 
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heat  in  a reactor  that  is  repeatedly  heated  and  cooled.  These  effects,  however,  are  ignored 
here,  resulting  in  a process  that  is  static  but  that  suffers  from  measurement  delay.  More 
specifically,  the  plant  open-loop  model  is  given  as 

yM(k)  = KPu(k-D)  (3-1) 

where  Kp  is  gain  of  the  plant,  D is  the  measurement  delay,  k is  an  indexing  notation  used 

to  number  batches  sequentially,  yM  (•)  is  the  measured  output,  and  u(-)  is  the  input. 

The  goal  of  this  study  is  to  provide  a comprehensive  theoretical  description  for 
the  single-input,  single  output  (SISO)  case  of  a static  run-to-run  plant  with  significant 
delay.  In  Section  3.2,  a suitable  proportional-integral  (PI)  controller  form  is  derived,  and 
compared  with  a standard  discrete  PI  controller.  Various  strategies  are  outlined  in 
Section  3.3  for  characterizing  the  parameters  of  a plant.  This  knowledge  is  required  for 
safe  and  effective  tuning  of  the  controller  for  closed-loop  feedback  control.  The 
application  of  a Smith  Predictor  is  investigated  in  Section  3.4  as  a possible  technique  to 
counter  the  degrading  effect  of  large  delays  on  closed-loop  performance.  It  is  shown  that 
even  linear  static  run-to-run  systems  with  significant  delay  are  deceptive  in  their 
simplicity,  as  a Smith  Predictor  can  lead  to  closed-loop  instability  for  very  small 
uncertainties  in  the  plant.  Finally,  Section  3.5  provides  conclusions  and  suggestions  for 
the  development  of  future  work. 

3.2  Derivation  of  a PI  Run-to-Run  Controller 

The  controller  for  a run-to-run  system  can  consider  only  past  data.  This  is  due  to 
the  fact  that  a batch  must  be  finished  completely  before  the  property  of  interest  can  be 
measured.  In  other  words,  the  output  y(k)  is  first  available  for  use  in  calculating  u(k+l). 
This  is  in  marked  contrast  with  the  situation  encountered  in  standard  sample-and-hold 
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data  systems,  where  it  is  assumed  that  the  output  y(k)  is  obtained  via  measurement  at  a 
sampling  instant  and  utilized  to  calculate  an  appropriate  control  input  u(k)  at  the  same 
sampling  instant.  For  the  run-to-run  system  in  question,  one  may  employ  a proportional- 
integral  control  law  that  calculates  the  input  from  knowledge  of  the  measured  error  e(k) 
obtained  from  the  processed  batch  k via  the  time-domain  control  law 


u(k  + l)  = Kce(k)  + KcRl[e(0)  + £(l)+...+8(k)]  , R,  >0  (3-2) 

where  Kc  is  the  proportion  gain  of  the  controller,  and  R,  is  the  reset  rate  of  the 

controller.  Taking  the  z-transform  yields  the  relationship 


u(k)z  = Kcs(k)  + KcR,[z“k  +z-k+,+...+z~1  +z°]£(k) 


(3-3) 


Recognizing  that  V z ' = — — and  rearranging  terms  as  needed  gives 

n Z-  1 


u(k)  = Kc  — £(k)  + KcRj  — — s(k)  (3-4) 

z z- 1 

or 

u(k)  = Kcf-  + RI-^|s(k)  (3-5) 

Vz  z -\J 

It  should  be  noted  that  R,  must  always  be  positive  since  a negative  value  for  reset  rate  is 

physically  meaningless.  Equation  (3-5)  is  similar  to  a standard  discrete  PI  controller  in 
position  form 


u(k)  = Kc 


1 + R,At 


s(k) 


(3-6) 


V Z-1 J 

where  At  is  the  sampling  period  of  the  controller  (Seborg  et  al.,  1989).  A main 
difference,  however,  is  that  (3-5)  has  an  extra  unit  delay  on  the  right-hand  side.  The  final 
subtle  difference  is  that  absolute  time  is  no  longer  a reference;  alternatively,  the  number 
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of  batches  marks  in  an  indirect  fashion  the  passage  of  time.  The  differences  arise  because 
a discrete  controller  can,  in  general,  be  applied  to  a continuous  plant,  which  has  both  an 
input  and  an  output  at  any  instant  in  time,  while  a discrete  plant  only  has  an  input  and 
output  value  for  each  batch  that  is  processed.  Appendix  C provides  a brief  review  of  how 
a discrete  controller  (Seborg  et  al.,  1989)  uses  the  current  error  to  implement  the  current 
input — a task  that  is  impossible  for  the  run-to-run  controller  to  accomplish. 

3.3  Plant  Stability  under  Feedback  Control 
Figure  3-1  shows  a block  diagram  of  the  closed-loop  system  considered 


Figure  3-1 . Generalized  block  diagram  for  a feedback  loop 
with  measurement  delay. 

1 1 3 


where  C = Kc[  - + R{  — — is  the  controller,  P = Kp  is  the  plant  gain,  M = z D is  the 

z - 1 


measurement  delay,  y is  the  output,  and  ym  is  the  measured  output.  Basic  block  diagram 


algebra  shows  that 


y 

r 


CP 


(3-7) 

i+CPM 

Upon  substitution  of  C,  P,  and  M,  closed-loop  stability  depends  upon  the  location 


of  the  uncancelled  poles  of  the  expression 
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y 


r 


i + : 


(3-8) 


When  rearranged  algebraically,  the  denominator  of  (3-8)  is  given  by  the  polynomial 


Appendix  D shows  that  (3-8)  never  has  unstable  pole-zero  cancellations  when  unstable 
poles  are  present,  provided  that:  (i)  R,  > 0,  and  (ii)  the  additional  constraints  required  by 
the  Jury  Array  are  met,  as  given  in  Section  3.3.1.  Therefore,  all  roots  of  (3-9)  that  lie 
outside  of  the  open  unit  disk  on  the  complex  plane  are  in  fact  unstable  poles  of  (3-8). 

Also,  inspection  of  (3-9)  reveals  that  there  are  only  three  lumped  parameters  that 
are  relevant  for  closed-loop  stability,  given  that  Kc  and  Kp  always  appear  as  the  product 
KCKP . Therefore,  the  locations  of  the  roots  of  equation  (3-9)  dictate  whether  or  not  the 
system  is  stable  for  a given  triplet  of  parameters  ( KCKP , R, , D ). 

3.3.1  Nominal  Stability  Region 

It  is  necessary  and  sufficient  for  the  poles  of  (3-9)  to  lie  within  the  unit  circle  to 
guarantee  stability,  based  on  the  analysis  provided  in  Appendix  D.  The  Jury  Test 
(Ogunnaike  and  Ray,  1994)  provides  criteria  in  two  steps  that  establish  necessity  and 
sufficiency  for  the  roots  of  a polynomial  to  lie  within  the  unit  circle.  The  first  step  of  the 
Jury  Test  states  that  given  the  polynomial 


where  the  leading  term  a0  is  positive,  there  are  three  necessary  conditions  for  the  roots  of 


P(z)  = zD+2  -zD+I  +[KCKP(R,  +l)]z  — KCKP 


(3-9) 


(3-10) 


a polynomial  to  lie  within  the  unit  circle: 
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1-  |a„|  < a0 

2.  P(l)  = Xa(>0 

i=0 

n (3-11,12,13) 

3.  P(— 1 ) = ^^  ("  1 )"  a;  > 0,  for  n even,  or 

i=0 

n 

P(-l)  = (“l)"  ‘ aj  < 0 for  n odd 

i=0 

Applying  (3-1 1 )-(3-l 3)  to  the  system  (3-9)  under  consideration  yields 

1.  |kck  I < i 

I c Pi 

2.  P(l)=  KcKpR,  > 0 (3-14,15,16) 

3.  P(-l)  = 2-KcKp(R, +2)>0,  for  D even,  or 
P(-l)  = -2  - KCKP(R,  +2)  < 0 for  D odd 

The  implications  of  this  set  of  necessary  conditions  are  immediate.  First,  from  the 

knowledge  that  R,  > 0,  based  on  (3-15),  it  is  necessary  for  both  KCKP  > 0 and  R,  > 0. 

Consequently,  because  of  (3-14),  it  is  necessary  that  0 < KCKP  < 1 . Both  conclusions  are 

reinforced  by  (3-16)  for  the  case  where  D is  even,  since  R,  > 0 implies  that  KCKP  < 1 . 

Finally,  given  R,  > 0,  the  case  where  D is  odd  implies  that  KCKP  > -1 , which  is  subset 

of  the  constraint  just  established  that  implied  0<  KCKP  < 1 . In  effect,  (3-14)  and  (3-15) 

can  be  reduced  to  yield  the  constraints 


1.  0 < KCKP  < 1 

2.  R,  >0 

3.  P(-1)  = 2-KCKP(R,  +2)>0,  forD  even,  or 
P(— 1)  = -2  — KCKP(R,  + 2)  <0  for  Dodd 


(3-17,18,19) 


The  second  step  in  applying  the  Jury  Test  requires  construction  of  the  Jury  Array. 


It  is  a table  composed  elements  that  are  calculated  from  the  coefficients  of  the  polynomial 
(3-10),  and  has  the  following  form: 


85 


where 


Row  1 

an 

an-l 

•"  a. 

Row  2 

ao 

ai 

•••  a„-. 

Row  3 

bn-l 

bn-2 

•••  bo 

Row  4 

bo 

b1 

•••  b„-l 

Row  2n-5 

P3 

P2 

Pi  Po 

Row  2n-4 

Po 

Pi 

p2  p3 

Row  2n  - 3 

q2 

q. 

q0 

Row  2n-2 

q0 

qi 

q2 

a„ 

an  i t 

bk  = 

n 

n-l-k 

; k = 

0,1,  ”,n 

ao 

ak+l 

Ok  = 

P3 

P2-k. 

k = 

0,1,2 

Po 

Pk+1 

(3-20) 


(3-21) 


The  elements  calculated  within  the  array  are  then  used  to  verify  whether  the 


following  additional  conditions  for  the  roots  of  the  polynomial  (3-10)  to  lie  within  the 


unit  circle  are  met: 


Ibn-l|  > |b0 


(3-22) 


|q2|  > |q0| 

When  applied  to  a polynomial  with  real  coefficients,  equations  (3-1 1)-(3-13)  and  (3-22) 
comprise  a set  of  necessary  and  sufficient  conditions. 

Example 

Let  D = 2 and  consider  P(z)  = z4  - z3  -t-[KcKp(R,  + l)]z-  KCKP . Since  D is 


even,  the  conditions  (3-1 7)-(3- 1 9)  reduce  to 
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1.  0 < KCKP  < 1 

2.  R,  > 0 (3-23,24,25) 

3.  2 - KCKP(R,  +2)  > 0 

Then  the  Jury  Array  for  the  problem  is  as  follows: 


Row  1 
Row  2 

-kckp 

1 

KcKp(RI+l) 

-1 

Row  3 

K2.K2-! 

-K2k2(Ri+1)+1 

Row  4 

KcKp-KcKp(RI+l) 

0 

Row  5 

(K2r2-1)2* 

(KcKp-KcKp(RI+l))2 

(KjKj-l)* 

(-K2k2(Ri+1)+1) 

Row  6 

(KcKp-KcKp(RI+l))* 

-(-K2k2(Ri+1)+1) 

(K2,K2-1)* 

(-K2k2(Ri+1)+1) 

0 -1  1 
0 KcKp(Rj+l)  -KcKp 

0 KcKp-KcKp(Rj+l) 

-K^Kp(Rj+l)+l  KCKp_1 

(KcKp-KcKp(RI+l))* 

-(-K2k2(Ri+1)+1) 

(K2,K2-1)2* 

(KcKp-KcKp(RI+l))2 


where  the  symbol  denotes  multiplication.  Therefore,  the  three  original  constraints  (3- 
17,18,19)  are  supplemented  by  two  additional  constraints  defined  by  (3-22)  for  this 
particular  example 


(3-26) 


|k2.k2-i|>|kckp-kckp(r1+i)| 

|(K2k2p-1)2(KcKp-KcKp(Ri+1))2|>|(KcKp-KcKp(R,+1))*-(-K(.K2p(Ri+1)+1)| 

Figure  3-2  represents  graphically  the  results  of  the  analysis  described  above  when 
it  is  repeated  for  various  values  of  D.  The  curves  are  the  nominal  stability  conditions  for 
different  sets  of  the  three  parameters  ( KCKP , R, , D ).  The  x-axis  represents  the  values  of 
R,,  and  the  y-axis  represents  the  product  KCKP . Each  curve  defines  the  stability  limit 
for  a fixed  value  of  the  measurement  delay  D , with  the  topmost  curve  corresponding  to 
D = 0,  the  curve  immediately  below  it  corresponding  to  D = 1 , and  so  on  down  to 
D = 20.  If  for  a given  D , the  point  defined  by  the  ordered  pair  (Rj  ,KcKp)  lies  below 

the  line,  the  system  is  stable  for  the  parameters  (Rj,KcK P,D).  A point  on  the  line  is 


critically  stable,  and  above  the  line,  unstable.  All  stability  regions  are  also  framed  by  the 
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x-axis  and  y-axis,  meaning  that  they  all  lie  entirely  within  the  first  quadrant,  as  dictated 
by  (3-1 7)  and  (3-18). 


R 


i 


Figure  3-2.  Nominal  stability  limit  curves  under  PI  control 
for  various  delay  values. 

These  curves  of  Figure  3-2  are  generated  using  the  root  finding  function  of 
Matlab  for  computational  efficiency,  but  confirmed  with  representative  Jury  Table  plots. 
Also,  for  clarity,  recall  that  the  controller  introduces  a unit  delay;  the  controller  delay  is 
separate  from  the  measurement  delay  D. 

3.3.2  Robustness  Analysis 

A static  plant  with  measurement  delay  of  the  form  (3-1)  under  PI  control  has  two 
parameters  that  are  known  exactly — the  gain  Kc  and  the  reset  time  R,  of  the  controller, 
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since  they  are  set  by  the  operator.  In  contrast,  it  is  possible  that  the  plant  gain  Kp  and 
the  measurement  delay  D may  not  be  known  precisely,  or  may  change  with  time.  More 
specifically,  the  delay  D may  be  uncertain,  but  is  restricted  to  be  a whole  number.  As  a 
result,  if  D changes,  the  stability  limit  simply  migrates  between  the  appropriate  curves  of 
Figure  3-2.  The  gain  Kp  of  the  plant,  however,  is  a continuous  parameter  that  may  be 
affected  by  error,  and  is  modeled  by  the  multiplicative  uncertainty 


with 


and 


Kp  = Kp(l  + s) 


£ < £ < £ 


(3-27) 


(3-28) 


0>  £ > -1 


(3-29) 


£ > 0 


(3-30) 


where  Kp  is  the  exactly  known  nominal  plant  gain,  and  £ is  a scalar  that  quantifies  the 
extent  of  modeling  error.  Equation  (3-27)  defines  the  unknown  actual  plant  gain  Kp  as  a 
multiple  of  the  exactly  known  nominal  plant  gain  Kp . The  range  of  possible  variation  is 
restricted  by  the  dual  inequality  (3-28),  where  £ is  the  lower  bound  and  £ is  the  upper 
bound.  Taken  together,  (3-29)  and  (3-30)  ensure  that  Kp  is  in  fact  one  of  the  possible 
values  of  Kp  (i.e.,  e = 0).  Finally,  Equation  (3-29)  further  restricts  the  error  so  that  Kp 
and  Kp  have  the  same  sign;  this  is  necessary  to  guarantee  (3-17)  once  the  value  of  Kc  is 


set. 
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Let  f(R,,D)  represent  an  arbitrary  curve  in  Figure  3-3.  The  robust  stability 
criterion  becomes  KCKP  < f(R,,D)  for  all  values  of  £ satisfying  (3-28).  Using  (3-27), 
this  condition  can  be  written  as 

KcKp(l  + s)  < f(R,,D)  (3-31) 

or  equivalently,  since  (1  + £)  > 0 

KrK°  < — — f(R, , D)  (3-32) 

c P (1  + e) 

Using  the  upper  constraint  bound  £ = £ yields  the  robust  stability  condition 

KrKp  < — l- — f(R,,D)  (3-33) 

c P (1  + e) 

In  other  words,  the  product  of  Kc  and  Kp  must  remain  below  the  upper  bound  given  by 
f(R,,D)  to  ensure  robust  stability.  The  nominal  stability  limit  is  obtained  for  the  case 
where  Kp  = K°  (i.e.,  £ = 0);  hence,  (3-33)  yields  the  nominal  stability  condition 

KcKp  < f(R,,D)  (3-34) 

Figure  3-3  shows  a comparison  of  the  boundaries  of  the  robust  stability  region 

given  by  (3-33)  and  the  nominal  stability  region  given  by  (3-9)  for  the  case  D = 5 and 

£ = 0.25 . Note  that  the  robust  stability  region  is  always  a constant  fraction  of  the 

nominal  stability  region:  namely,  the  fraction  — that  results  from  taking  the  ratio  of 

(l  + £) 

the  right-hand  sides  of  (3-33)  and  (3-34). 

The  robust  stability  limit  is,  in  effect,  a gain  margin  that  reduces  the  ordinate  of 
each  critically  stable  ordered  pair  by  the  factor  1/(1  + e) . The  graph  resulting  from  a 
value  of  e = 1 / 9 is  shown  in  Figure  3-4.  Note  that  the  graph  has  the  same  general  shape 
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Figure  3-3.  Representation  of  the  boundary  robust  stability 
region  given  by  (3-33)  (solid  line)  and  the  boundary  of  the 
nominal  stability  region  given  by  (3-34)  (dashed  line)  for  a 
delay  D = 5 . 

as  that  shown  in  Figure  3-2  for  the  nominal  stability  region,  except  that  each  y-value  has 

been  scaled  by  the  factor  ^ = 0.9. 

l + (l/9) 

Given  a static  plant  of  the  form  (3-1)  including  a measurement  delay  D,  a key 
engineering  task  is  to  identify  the  error  bounds  in  the  modeling  parameters  Kp  and  D.  A 

systematic  design  of  experiments  (Box  el  al.,  1978)  may  be  an  effective  strategy  for 
determining  the  range  of  gain  values  Kp , or  equivalently,  finding  bounds  on  the  error 
parameter  s.  In  fact,  standard  linear  regression  studies  (Kreyszig,  1999)  using 
experimental  data  can  be  used  to  yield  confidence  intervals  directly  related  to  the 
uncertainty  parameter  £ . In  addition,  if  the  process  is  already  in  operation,  historical  data 
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Figure  3-4.  Robust  stability  limit  curves  under  PI  control 
for  e = (1/9)  and  various  delay  values. 

may  be  used  to  help  reverse-engineer  the  range  of  the  plant  gain  values. 

Identifying  the  value  of  D is  generally  more  straightforward.  Whether  in  the 
design  stage  or  already  in  operation,  it  should  be  trivial  to  count  the  number  of  units  that 
are  en  route  from  the  point  at  which  the  input  value  is  applied  to  the  point  at  which  the 
output  is  measured.  However,  it  should  be  noted  that  it  is  possible  for  the  delay  to  vary 
since  not  all  unit  operations  are  necessarily  occupied  by  a batch  at  a given  time.  The 
importance  of  this  scenario  is  discussed  in  detail  in  Section  3.4. 

3.3.3  Tuning  Correlations 

Generating  the  curves  represented  by  f(R,;D)  is  computationally  expensive, 
since  the  polynomials  generated  by  the  Jury  Table  are  of  order  (D+l)  . It  is  observed  that 
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they  resemble  the  frequency  response  for  a sinusoidal  input  (Levine,  1 996),  and  can  be 
approximated  by  a two-parameter  (a,p)  model 


f'(R„D) 


1 

(aR,)4  + 2(aR,)2(2p- 1)  + 1 


(3-35) 


Hence,  f'(R,,D)  = f(R,,D)for  appropriate  values  of  parameters  a and  P . As  a 
model  with  only  two  parameters,  the  optimal  values  of  a = a(D)  and  p = P(D)  can  be 
found  easily  through  a computationally  inexpensive  exhaustive  numerical  search.  This  fit 
can  be  optimized  at  the  discretion  of  the  user  with  least  squares,  minimum  absolute  error, 
etc.  Results  are  given  in  Appendix  E for  the  case  of  minimum  absolute  error,  and  the 


Figure  3-5.  Optimal-fit  parameters  for  the  model  (3-35) 
using  the  minimum  absolute  error  criterion:  (a)  parameter 
a , (b)  parameter  p . 
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resulting  values  of  a and  p as  a function  of  D are  given  in  Figure  3-5.  In  addition, 
Figure  3-6  shows  the  plots  of  the  exact  curves  KcKp  = f(R,,D)  and  of  the  approximate 
curves  KcKp  = f(R,  ,D) . Note  that  only  the  curves  for  delay  values  of  2,  5,  8,  11,  14,  and 

1 7 are  plotted;  this  shows  the  performance  of  the  fit  over  a wide  range  of  values  of  D, 
while  retaining  graphical  clarity. 


Figure  3-6.  Visual  comparison  of  actual  stability  curves 
with  those  generated  by  2-parameter  correlation. 

Figure  3-6  shows  that  some  of  the  approximations  do  not  provide  a close  fit,  but 
the  approximate  curves  are  nonetheless  useful  when  considering  how  tuning  heuristics 
are  implemented.  For  instance,  a common  method  is  to  begin  along  the  y-axis  and  select 
one  half  of  the  maximum  value  allowed,  in  this  case  KCKP  = 0.5  for  all  delays.  The 
abscissa  of  the  point  found  in  this  fashion  is  the  recommended  value  of  KCKP.  Next, 
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select  the  relevant  delay  curve  for  the  problem,  and  move  out  one  third  of  the  way  to  the 
curve  in  the  x-direction  while  keeping  KCKP  = 0.5  . This  resulting  point  is  the  suggested 
initial  tuning.  With  such  gross  approximations,  it  should  be  clear  that  the  small  amount 
of  error  introduced  using  the  two-parameter  fit  is  negligible.  Moreover,  with  suggested 
tunings  so  far  from  the  limits  of  stability,  robustness  is  typically  ensured. 

3.4  Plant  Stability  under  PI  Control  with  a Smith  Predictor 
The  application  of  a Smith  Predictor  (Marlin,  2000)  for  deadtime  compensation  is 
a potential  remedy  for  a system  affected  with  significant  delay  values.  When  added  to 
the  system  discussed  in  Section  3.3,  the  closed-loop  diagram  becomes  as  illustrated  in 
Figure  3-7 


Figure  3-7.  Generalized  block  diagram  for  a feedback  loop 
of  a static  plant  and  measurement  delay  using  a PI 
controller  with  a Smith  Predictor. 

As  discussed  in  Section  3.3,  the  plant  P and  the  measurement  system  M have  two 
uncertain  parameters:  namely  Kp  and  D.  To  implement  the  Smith  Predictor,  estimated 

values  Kp  and  D°  are  required  to  build  the  Smith  Predictor  operators 

P”  = Kp  and  M°  = z‘D°.  It  is  desirable  to  investigate  the  robustness  of  the  Smith 

Predictor  scheme  with  respect  to  errors  in  the  estimates  of  the  parameter  Kp  and  D. 
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3.4.1  Closed-loop  dynamics  under  PI  control  with  a Smith  Predictor 

Standard  block  diagram  algebra  on  Figure  3-7  yields  the  transfer  function  for  the 
inner  loop 


e 1 + CP°  -CM°P° 


(3-36) 


Using  (3-36)  in  the  relationship  y = Pu  and  recognizing  that  e = r - y yields 


y = P 


C 


OnO 


U + CP  -CMUP 

which  can  be  readily  rearranged  into  the  form 


(r  - My) 


(3-37) 


y CP 

r ~ 1 + CP°  - CM°P°  + CMP 


(3-38) 


As  in  the  case  considered  in  Section  3.3,  the  stability  region  is  determined  by  the  roots  of 
the  characteristic  polynomial  P(z)  = 1 + CP0  -CM°P°  +CMP,  as  shown  in  Appendix  F. 
Substituting  the  values  for  C,  M,  and  P defined  in  Section  3.3,  as  well  as  those  for 
P°  and  M1'  defined  in  Section  3.4  yields  the  characteristic  equation 


1 + K, 


-+R, 

k z 


z — 


Kp  - Kc 


f 

<z 


- + R1-^-|k0dz“d" 


z — 


+Kr 


-+R, 

Vz 


z-  ] 


Kpz“  = 0 


(3-39) 


Multiplying  each  term  by  z(z  - l)zDzD  and  grouping  like  terms  yields 


zD"+D+2  +[KcKp(R,  + 1) - l]zD"+D+l  - (KcKp)zD"+D 
-[KcKp(R,  + l)]zD+1  +(KrK°p)zD 
+[KCKP(R,  + l)]zD  +1  - (KcKp)zD  = 0 


(3-40) 
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Provided  P = P°  and  M = M°,  inspection  reveals  that  the  last  four  terms  of  (3-40)  sum  to 
zero,  leaving  the  simplified  nominal  characteristic  equation 
z2  +[KcKp(R,  +l)-l]z-KcKp  =0,  which  is  equivalent  to  (3-9)  for  the  special  case  of 
an  undelayed  plant  (i.e.,  D=0). 

3.4.2  Stability  Region  as  a Function  of  Model  Uncertainty 

As  the  true  plant  parameters  are  allowed  to  deviate  from  their  estimations,  the 
loop  may  become  unstable.  For  illustration  purposes,  the  stability  space  is  depicted  in 
Figure  3-8  for  the  set  of  conditions  KCKP  = 0.5,  R,  = 0.2,  and  D°  = 6.  At  zero  on  either 
axis,  the  estimated  parameter  matches  the  true  parameter  exactly.  Moving  away  from 
zero  represents  increasing  the  mismatch,  with  negative  values  representing 
underestimates  and  positive  values  representing  overestimates. 

For  example,  in  the  case  where  D°  - D = 0 , the  stability  range  for  KCKP  is 
approximately  -0.45  < KC(KP  - Kp)  < 0.48  . In  contrast,  at  D°  - D = ± 1 , there  are  no 
stabilizing  KCKP  values.  Similarly,  for  D°  - D = + 2 , the  system  remains  stable  for  the 
approximate  range  0.1  < KC(KP  - Kp)  < 0.45,  and  so  forth  for  greater  values  of 
mismatch.  Note  that  the  space  is  not  symmetrical,  which  is  most  dramatically  shown  by 
a significant  stability  range  for  D°  - D = 5 , where  none  exists  for  D°  - D = -5 . In 
addition,  for  the  subset  of  cases  in  which  Kp  = Kp , the  only  stable  value  of  D°  -D  is 
zero.  Therefore,  for  all  practical  purposes,  it  is  essential  to  know  the  value  of  the  delay 
precisely  in  order  to  consider  using  a Smith  Predictor  on  this  type  of  system. 

To  interpret  this  diagram,  first  recognize  that,  while  not  explicitly  shown,  there 
are  physical  limits  to  underestimates  (negative  numbers),  and  no  bound  to  overestimates 
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D°-D 


Figure  3-8.  Stability  region  for  a run-to-run  controller 
applied  to  a discrete  plant.  The  exact  nominal  parameters 

are  KcKp  = 0.5,  R,  = 0.2,  and  D°  = 6 and  uncertain 
parameters,  plant  gain,  Kp , and  measurement  delay,  D , 
are  allowed  to  vary. 

(positive  numbers).  Second,  the  delays  must  be  whole  numbers  (and  therefore  errors  in 
them  must  be  whole  numbers),  while  errors  in  gains  are  differences  of  real  numbers. 
That  is  why  the  stability  space  only  shows  discrete  values  across  the  x-axis.  In  addition, 
the  stability  space  is  continuous  vertically  within  each  grouping  of  circles;  they  are 
spread  out  only  for  clarity.  The  most  important  feature  of  the  diagram  is  the  intolerance 
of  the  system  to  error  of  as  little  as  one  unit  in  delay.  Furthermore,  at  any  error  in  delay 
of  other  than  zero,  the  system  in  unstable  for  a perfect  match  of  the  gain  and  its  estimate. 
Clearly,  unless  the  engineer  has  complete  confidence  in  the  value  and  constancy  of  the 
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delay,  a Smith  Predictor  will  result  in  unstable  behavior  even  for  little  or  no  error  in  the  in 
the  gain  estimate. 

For  contrast  purposes,  Figure  3-9  shows  the  analogous  stability  space  for  a 
discrete  controller  applied  to  a continuous  system.  Specifically,  the  pertinent  simulation 

parameters  are  KcKp  = 0.5,  R,  = 0.2,  and  D°  = 6 . The  area  is  represented  by  groups  of 
crosses,  and  is  continuous  vertically;  they  are  spaced  only  for  clarity.  The  circles 
represent  the  same  stability  space  depicted  in  Figure  3-8.  It  is  clear  that  the  stability 
region  for  the  discrete  controller  applied  to  a continuous  system  is  much  larger  and  not 
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Figure  3-9.  Stability  region  for  a discrete  controller  applied 
to  a continuous  plant.  The  exact  nominal  parameters  are 
KcKp  = 0.5,  R,  = 0.2,  and  D°  = 6 and  uncertain 
parameters,  plant  gain,  Kp , and  measurement  delay,  D , 
are  allowed  to  vary. 


99 


fragmented,  unlike  the  similar  run-to-run  controller/discrete  plant  stability  region  ( i.e .,  the 
circles).  Most  notably,  it  shows  that  a continuous  system  with  a discrete  controller  can 
tolerate  errors  in  the  estimation  of  the  delay  for  every  value  within  the  range  from 
D°  - D = -6  to  D°  - D = 6 when  there  is  no  error  present  in  the  plant  gain  value. 
Furthermore,  the  discrete  controller  system  is  stable  for  the  approximate  range 
-1.4  < Kc(Kp  - Kp)  < 0.4  when  the  error  in  delay  is  zero.  Recall  that  that  is  not  the  case 
for  the  run-to-run  controller/discrete  plant  system.  In  other  words,  the  discrete  controller 
with  Smith  Predictor/continuous  plant  system  depicted  in  Figure  3-9  is  robust  enough  to 
tolerate  small  errors  in  either  of  the  two  variable  parameters,  while  the  run-to-run 
controller  with  Smith  Predictor/discrete  plant  system  is  not. 

In  summary,  while  the  two  types  of  systems  may  appear  very  similar,  there  are 
fundamental  differences  that  significantly  affect  the  performance  and  stability  space  of 
each  configuration.  Therefore,  it  is  important  that  a control  designer  recognize  the 
distinction  between  the  two  subtly  difference  arrangements.  If  the  system  is  continuous, 
a Smith  Predictor  is  a potential  candidate  to  remedy  the  adverse  effects  of  measurement 
delay.  However,  it  is  concluded  from  the  robustness  analysis  (Section  3.3.2)  and  the 
tuning  correlations  (Section  3.3.3)  that  PI  control  is  a safer  strategy  for  a run-to-run 
controller  applied  to  a discrete  system,  despite  the  fact  that  limited  performance  can  be 
expected  because  of  the  significant  delay. 

3.5  Conclusions 

The  PI  form  of  the  run-to-run  controller  has  been  derived.  In  addition,  tuning 
heuristics  have  been  provided  for  a wide  range  of  possible  plant  gain  and  delay  values. 
With  them,  stable  plant  behavior  can  be  guaranteed,  although  limited  performance  must 
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be  tolerated.  However,  it  has  been  shown  that  stability  cannot  be  guaranteed  when 
applying  Smith  Predictors  to  run-to-run  static  systems  with  measurement  delay  in  the 
presence  of  even  very  small  modeling  error.  The  implication  of  Figure  3-8  is  that  the 
root  locus  of  the  characteristic  equation  for  the  system  with  a Smith  Predictor  moves 
around  as  a function  of  D°-D  and  Kc(Kp-Kp).  Furthermore,  as  D changes,  it  is 
apparent  that  one  or  more  roots  migrate  into  and  out  of  the  unit  circle  (the  critical  limit  of 
stability).  Therefore,  it  is  recommended  that  simple  PI  control  be  used  for  run-to-run 
static  systems  with  significant  delay.  Furthermore,  a set  of  tuning  heuristics  are  provided 
that  will  yield  satisfactory  performance. 

A comprehensive  study  of  the  behavior  of  the  roots  of  the  characteristic 
polynomial  may  be  helpful  in  finding  ways  to  modify  the  controller  to  be  more  robust 
with  a Smith  Predictor.  Without  that  development,  it  is  necessary  to  use  more  traditional 
control  approaches,  even  though  performance  will  be  limited.  A possible  starting  point 
may  be  to  study  the  discrete  analog  of  the  findings  of  (Adam  et  al.,  2000),  who  have 
shown  that  discrete  stability  ranges  exist  for  a Smith  Predictor  applied  to  a continuous 
system  with  an  uncertain  time-delay  parameter.  In  their  work,  the  frequency  response  in 
the  Nyquist  plane  is  shown  to  rotate  about  the  origin  as  the  delay  value  varies,  which  can 
alter  the  number  of  encirclements  of  the  critical  point  (Latchman  and  Crisalle,  1995), 
which  is  a parameter  that  is  central  to  the  stability  of  a continuous  system. 


CHAPTER  4 

CONTROL  OF  pH  USING  A NOVEL  OPTICAL  SENSOR 

4.1  Introduction 

The  pH  of  a solution  can  have  extreme  effects  on  the  environment.  For  example, 
river  ecology  is  extremely  sensitive  to  pH,  and  manufacturing  facilities  that  discharge 
materials  into  river  water  must  be  careful  to  not  create  a toxic  environment.  In  many 
chemical  reactions,  the  pH  determines  both  extent  of  reaction  and  selectivity,  and  hence  it 
can  be  a determining  factor  for  the  profitability  of  a design  (Fogler,  1992).  Finally,  pH  is 
also  an  important  attribute  of  the  final  product  for  goods  ranging  from  car  batteries  to  soft 
drinks. 

Measurement  of  pH  is  a commonly  accepted  exercise,  but  the  true  range  of  the 
measurement  demands  involved  is  unusually  large.  The  definition  of  pH  is 

pH  = -log[H  + ] (4-1) 

where  [H+]  is  the  hydrogen  ion  concentration  in  moles  per  liter  (mol/L  or  M).  More 

accurately,  it  is  a function  of  hydrogen  ion  activity,  but  the  two  are  essentially  equivalent 

for  the  purposes  of  pH  measurement  in  most  systems  (McMillan,  1 994).  Furthermore,  it 

is  believed  that  hydrogen  ions  exist  as  hydronium  ions  [H30+]  in  aqueous  solution;  this 

can  also  be  ignored  for  simplicity.  Therefore,  the  problem  of  controlling  pH  is  really  the 

problem  of  controlling  concentration  of  the  chemical  species  it  represents,  namely, 

hydrogen  ion  concentration.  Even  when  considering  an  operating  range  of  2 < pH  < 12 

(more  extreme  values  are  not  uncommon),  this  constitutes  a tremendous  range  (0.01  to 
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0.000000000001  M)  of  hydrogen  ion  concentration  conditions  that  must  be  taken  into 
account.  Additionally,  some  electrodes  can  respond  to  changes  of  as  little  as  .001  pH 
units,  and  under  certain  conditions  at  neutrality  (pH  = 7) , that  corresponds  to  a hydrogen 
ion  concentration  change  of  less  than  0.000000001  M (McMillan,  1994).  The 
combination  of  range  and  sensitivity  demands  makes  it  a major  challenge  to  design  a 
control  scheme  that  that  can  work  for  such  varied  conditions  with  acceptable  precision 
due  to  the  fundamental  laws  governing  signal-to-noise  ratio  (Williams,  1991). 

From  the  viewpoint  of  control  engineering,  the  plant  (a  solution  with  an  arbitrary 
pH  value  residing  in  a vessel  or  flowing  in  a stream)  is  highly  non-linear  over  its 
operating  range.  The  chemistry  involved  in  acid/base  reactions  is  a major  underlying 
reason  for  this  trait  (see  Appendix  G).  Furthermore,  when  considering  a system 
including  at  least  one  weak  acid  or  base,  the  amount  of  additional  acid  or  base  needed  to 
bring  a volume  of  fluid  to  a desired  pH  varies  according  to  the  strengths  of  the  acid(s)  or 
base(s),  their  dissociation  constants,  and  the  initial  pH.  In  addition,  pH  regions  exist 
where  the  fluid  volume  resists  a change  in  pH  when  additional  acid  or  base  is  added — a 
phenomenon  referred  to  as  buffering  (Zumdahl,  1989).  In  other  words,  different  fluids  of 
the  same  volume  and  pH  can  require  significantly  different  amounts  of  acid  or  base  to 
achieve  the  same  shift  in  pH,  and  proceed  along  the  transition  with  very  different 
trajectories.  Buffering  contributes  significantly  to  an  underlying  difficulty  of  pH  control, 
namely,  that  the  gain  of  a pH  system  can  vary  by  several  orders  of  magnitude  from  point 
to  point  throughout  the  operating  range. 

Finally,  in  most  realistic  pH  control  applications,  the  acidic  and  basic  components 
of  the  stream  that  must  have  its  pH  adjusted  changes  over  time.  This  introduces  different 
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dissociation  constants,  and  therefore  makes  the  already  highly  variable  gain  (with  respect 
to  operating  range)  time-varying  as  well.  Even  though  this  phenomenon  can  still  be 
controlled  in  the  traditional  manner  (feedback  control  strives  to  eliminate  offset),  the 
continuous  variations  often  severely  degrade  control  performance,  especially  when  the 
tuning  parameters  are  static.  The  complicated  nature  of  describing  the  pH  of  a 
combination  of  streams  has  inspired  considerable  work  in  the  area  of  open-loop  pH 
reactor  modeling  (McAvoy,  1972;  McAvoy  et  al.,  1972;  Gustafsson  et  al.,  1995). 

There  have  been  many  attempts  to  devise  schemes  for  closed-loop  control  that  are 
effective  despite  the  problems  of  range,  sensitivity,  and  plant  gain.  These  include  the 
implementation  of  dual  base  injectors  (Williams  et  al.,  1990),  and  the  use  of  a smaller, 
separate  reactor  (Sung  and  Lee,  1 995)  to  quickly  identify  the  behavior  of  the  system  and 
use  that  knowledge  to  improve  control  performance.  Furthermore,  advanced  control 
strategies  such  as  adaptive  control  (Mahuli  et  al.,  1993;  Sung  et  al.,  1998)  have  been  used 
to  ensure  that  the  controller  remains  optimally  tuned  for  the  observed  operating 
conditions.  Finally,  it  has  been  proposed  that  using  alternative  control  variables  may  be 
an  effective  way  to  improve  pH  control  (Gustafsson  and  Waller,  1983;  Waller  and 
Makila,  1981;  Narayanan  et  al.,  1998).  In  one  way  or  another,  these  arrangements 
attempt  to  produce  a system  that  is  easier  to  control  despite  the  problems  stemming  from 
range,  sensitivity,  and  variable  gain. 

One  common  feature  of  most  methods  is  that  they  use  a standard  electrode  as  a 
measurement  device  to  close  the  control  loop.  Because  of  the  physics  that  govern 
electrodes  (see  Section  4.2.2),  the  electronic  signal  sent  to  the  controller  is  proportional  to 
pH.  This  introduces  an  additional  and  significant  non-linearity  to  the  system,  since  [H+] 
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is  the  chemical  species  being  controlled.  An  optical  pH  sensor  (Janowiak,  2003), 
commonly  referred  to  as  an  optrode  (see  Section  4.2.3)  being  developed  at  the  University 
of  South  Florida  generates  a signal  proportional  to  [H+],  not  pH,  and  therefore  holds 
potential  for  eliminating  the  non-linearity  of  the  feedback  portion  of  the  loop.  While  not 
addressing  the  non-linear  physics  that  govern  the  chemistry  in  the  plant,  the  optical 
sensor  with  a signal  proportional  to  the  true  variable  of  interest  ([H+])  could  affect  the 
closed-loop  dynamic  behavior  of  the  system  dramatically.  The  goal  of  this  study  is  to 
determine  if  a linear  optical  pH  sensor  achieves  feedback  control  performance  superior  to 
that  of  the  traditional,  highly  non-linear,  electrode. 

There  are  several  components  required  to  achieve  the  project  goal.  First,  an 
accurate  model  must  be  developed  to  simulate  a pH  control  system;  Section  4.2  discusses 
the  process  in  detail.  Section  4.3  outlines  a simulation  study  designed  to  access  the 
relative  performance  of  the  two  probes  over  a range  of  response  characteristics,  both 
without  noise  and  with  realistic  noise  injected.  Section  4.4  addresses  issues  pertinent  to 
the  construction  of  a pilot-scale  pH  control  system  for  the  purpose  of  experimental  trials. 
Conclusions  are  drawn  about  the  relative  effectiveness  of  the  two  sensors  with  respect  to 
dynamic  control  in  Section  4.5. 

4.2  Modeling  and  Simulation  of  a pH  Control  System 

There  are  three  primary  components  to  a basic  feedback  control  scheme.  Their 
relative  arrangement  is  represented  in  the  abstract  block  diagram  that  is  Figure  4-1,  where 
C represents  the  controller,  P is  the  plant,  M is  the  measurement  device,  r is  the  system 
setpoint,  u is  the  input  to  the  system,  and  y and  ym  are  the  output  and  measured  output  of 
the  system,  respectively.  Each  of  the  three  components  C,  P,  and  M is  modeled 
separately  before  the  dynamic  response  of  the  closed  system  is  evaluated. 
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Figure  4-1.  Abstract  block  diagram  of  a simple  feedback 
control  system. 

When  modeling  the  plant  P,  it  is  noted  that  chemical  reactions  among  acids  and 
bases  are  relatively  fast  (Zumdahl,  1989).  Therefore,  using  a continuously  stirred  tank 
reactor  (CSTR)  is  a plausible  equipment  choice  for  pH  control.  It  is  assumed  that 
sufficient  agitation  allows  any  reaction  to  proceed  to  equilibrium  well  before  the  fluid 
exits  the  reactor.  Additionally,  the  simplicity  of  the  apparatus  helps  minimize  the  number 
of  parameters  associated  with  the  entire  control  system.  As  a result,  a more  confident 
determination  can  be  made  of  the  relative  advantages  and  disadvantages  of  using  an 
optical  pH  sensor  in  lieu  of  a traditional  electrode.  Section  4.2.1  is  a step-by-step 
description  of  the  modeling  process  adopted  to  describe  a CSTR  for  pH  control. 

The  most  common  control  instrument  C used  for  pH  control  systems  is  the 
proportional-integral  (PI)  controller 


u(s)  = Kc 


1+  — 


e(s) 


(4-2) 


v w 

where  Kc  is  the  gain  of  the  controller,  x,  is  the  integral  time  constant  of  the  controller,  e 


is  the  error,  and  u is  the  input  to  the  CSTR.  Derivative  action  is  not  used,  since 
measurement  noise  typically  degrades  the  response  to  a point  where  the  benefits  of 
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derivative  action  are  significantly  outweighed  by  the  deleterious  effects  of  noise.  In 
addition,  more  advanced  control  schemes,  such  as  cascade  or  Smith  Predictor 
configurations,  are  assumed  to  be  unnecessary.  They  are  generally  used  in  situations 
where  input  disturbances  or  long  transport  delays  degrade  the  quality  of  control 
performance,  problems  not  associated  the  chosen  CSTR  model. 

The  control  loop  is  closed  with  the  addition  of  a measurement  device  M.  Both  the 
electrode  and  optrode  work  by  continuously  undergoing  chemical  reaction;  see  Sections 
4.2.2  and  4.2.3  for  a detailed  discussion  as  well  as  a mathematical  model  for  each  device. 
Because  the  physics  of  the  two  instruments  are  different,  the  closed  loop  response  can  be 
expected  to  be  different  as  well.  Establishing  a quantitative  criterion  to  evaluate 
performance  should  indicate  whether  or  not  there  is  a significant  difference  in  the 
behavior  of  the  closed-loop  system  when  utilizing  one  sensor  versus  the  other. 

4.2.1  Process  Model  of  pH  Control  in  a CSTR 

The  process  model  consists  of  a mixing  tank  and  three  streams — two  entering  the 
CSTR  and  one  leaving  it,  as  shown  in  Figure  4-2.  One  of  the  inbound  streams  is  the 
influent,  and  typically  consists  of  a mixture  of  a number  of  weak  acids  that  is  typical  of 
industrial  process  by-products.  The  other  incoming  stream  is  the  reagent,  which  is  a 
strong  base  used  to  neutralize  the  acids.  The  effluent  should  have  the  desired  pH  for  the 
next  destination  of  the  stream.  Note  that  in  a production  process,  the  influent  may  change 
from  acidic  to  basic,  or  back  again.  A control  system  for  such  a process  would  require  an 
additional  input  stream  of  acidic  reagent.  Only  the  case  of  an  acidic  stream  being  treated 
with  a basic  reagent  is  considered  in  this  study  because  this  scenario  is  no  more 


illustrative  in  sensor  comparison. 
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A schematic  of  the  reactor  with  volume  V used  in  this  study  is  shown  in  Figure  4- 
2.  The  influent  stream  with  flowrate  Q,  is  chosen  to  be  acetic  acid  of  concentration 
[A] jn,  a commonly  available  weak  acid  with  a well-known  dissociation  constant.  Second, 
the  reagent  stream  with  flowrate  Q2  is  chosen  to  be  sodium  hydroxide  (also  commonly 
available)  of  concentration  [NaOH]jn.  Note  that  the  following  inlet  concentration  holds: 
[NaOH]jn  = [Na+]in=  [OH~]jn.  Finally,  the  effluent  stream  leaving  the  tank  at  a flowrate 
Q3  is  a mixture  of  both  acetic  acid  and  sodium  hydroxide  in  concentrations,  [A]  and 
[Na+],  respectively.  The  flowrates  Q,,Q2,  and  Q,  are  measured  in  liters  per  second 
(L/s),  and  [ ] denotes  a concentration  of  species  ( ) in  moles  per  liter  (mol/L  or  M). 


Figure  4-2.  Schematic  of  a CSTR  used  for  pFI  control. 

The  entire  system  is  described  by  a set  of  six  equalities,  of  which  two  are 

differential  equations  and  four  are  algebraic  equations 


d[A] 


Qi  + Q2 

V 


[A]  + (A  [A]„ 


dt 


(4-3) 
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d[Na+] 

dt 


= - "Q|  + Q-2~l[Na+  ] + { — l[Na+  ]in 

l v / lv/ 

(4-4) 

Y [H+][A“] 
3 [HA] 

(4-5) 

Kw=[H+][OH-] 

(4-6) 

[A]  = [HA]  + [A-] 

(4-7) 

[Na+]  + [H  + ] = [A  ] + [OH  ] 

(4-8) 

An  important  feature  of  the  system  is  that  the  effluent  stream  Q3  is  created  by 
spillover  from  the  CSTR,  ensuring  that  the  volume  V is  constant.  Assuming  that  all  fluid 
streams  are  aqueous  with  approximately  the  same  density,  a total  mass  balance  yields 


”~ApQi  +pQ;  -pQ,  (4-9) 

dt 

Assuming  that  the  density  p is  constant  and  using  the  fact  that  V is  also  constant,  it 
follows  that 


Q,+Q2=Q3  (4-10) 

The  two  differential  equations  represent  chemical  species  that  do  not  change  in 
the  reactor,  and  are  derived  via  mass  balances.  The  first  is  total  acetic  acid  concentration, 
denoted  [A],  which  is  the  sum  of  the  concentrations  of  the  protonated  and  unprotonated 

molecules,  i. e.  [A]  = [HA]  + [ A ” ] . 


d([A]V) 

dt 


= [A]inQ,-[A]Q3 


(4-11) 


Now  using  the  fact  that  V is  constant  and  the  flowrate  equality  (4- 1 0) 
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V^l  = [A],„Q,-[A](Q,+Qi) 
dt 


(4-12) 


or 


d[A]  (Q]+Q 


dt 


V 


[A]  + 


f |[A]in 


which  is  the  first  of  the  six  equations  (4-3)  that  describe  the  system.  The  second  species 
balance  required  is  for  sodium  ions,  yielding 


d([Na+]V) 

dt 


= [Na+]inQ2-[Na+]Q3 


which  can  be  manipulated  as  before  to  give 


(4-13) 


d[Na+]_  f Q,  +Q2 


\ f 

+ 


[Na’I  + ^SljlNa*],, 


dt  V V 

or  the  second  of  the  six  equations  (4-4)  that  describe  the  system. 

Equations  (4-3)  and  (4-4)  define  the  two  states  of  the  reactor.  Because  there  are  a 
total  of  six  unknowns 


[A],[Na+],[H+],[A^],[HA],  and  [OET] 

four  additional  equations  are  required  to  account  for  all  of  the  relevant  species.  First, 
both  acetic  acid  and  water  have  dissociation  constants,  Ka  and  Kw,  which  are  reported  as 
(4-5)  and  (4-6)  respectively,  such  that 


and 


Ka:  = 


[H+][A~] 

[HA] 


so  that 


Kw:=[H  + ][OH-] 
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[OH~]  = -^-  (4-14) 

[H  + ] 

Second,  there  is  the  aforementioned  total  acid  concentration  (4-7) 


[A]  = [HA]  + [A“] 

Third,  charge  neutrality  is  required  for  a mixture  of  neutral  liquids,  and  is  given  by  (4-8) 
in  the  form 


[Na+]  + [H  + ] = [A~]  + [OH~] 

The  quantity  of  interest,  hydrogen  ion  concentration,  is  found  through  repeated 
substitutions.  Rearranging  (4-5)  in  the  form 


KJHA]  = [H  + ][A  ] (4-15) 

and  then  rearranging  (4-7)  to  find  that 

[HA]  = [A]-[A-]  (4-16) 

yields  after  (4-16)  is  substituted  into  (4-15)  and  a rearrangement  of  terms 


Ka[A]  = [H  + ][A-]  + Ka[A-] 


Solving  now  for  the  concentration  of  the  anion 


(4-17) 


K„[A] 
[H+]  + Ka 


(4-18) 


Now  using  (4-18)  and  (4-14)  in  the  charge  neutrality  equation  (4-8)  and  rearranging 


terms  yields 


[Na+]  + [H  + ]- 


Ka[A] 
[H+]  + Ka 


*-  = 0 
[H+] 


(4-19) 


or 
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[Na+ ]([H  + ] + Ka  )[H+ ] + ([H + ] + Ka  )[H+ ]2 
-Ka[A][H  + ]-([H  + ] + Ka)Kw=0 

which  is  easily  cast  in  the  form  of  a third-order  polynomial 


[H  + ]3  +([Na+]  + Ka)[H  + ]2  + 
(Ka[Na+]-Ka[A]-Kw)[H  + ]-KaKw  =0 


(4-21) 


It  is  critical  to  select  the  proper  root  of  the  polynomial  (4-21)  to  find  [H+],  Using 
Descartes'  Rule  of  Signs  (Borowski  and  Borwein,  1991)  to  analyze  the  coefficients 
proves  that  Equation  (4-21)  has  exactly  two  negative  real  roots  and  one  positive  real  root 
(see  Appendix  H).  Because  negative  concentration  values  have  no  physical  meaning,  the 
correct  solution  is  given  by  the  positive  root. 

4.2.2  Traditional  Electrodes 
4.2.2. 1 Theory  of  operation 

The  basic  design  of  traditional  pH  electrode  sensors  consists  of  a coated  glass 
electrode  filled  with  a fluid  of  known  pH.  The  coated  glass  is  pH  sensitive,  and  if 
immersed  in  a fluid,  galvanic  half-cell  potentials  develop  across  the  external  fluid-glass 
boundary  and  the  glass-fill  solution  boundary  (McMillan,  1994).  If  care  is  taken  to 
ensure  that  the  hydrated  gel  layers  on  each  side  of  the  glass  are  suitably  similar,  and  that 
sufficient  free  ions  are  available  to  migrate,  the  two  potentials  can  be  compared  directly. 
The  resulting  voltage  difference  is  a function  of  both  pH  and  temperature. 

The  ability  of  a standard  glass  electrode  to  measure  pH  is  rooted  in  the  free 
energy  present  in  the  system,  and  is  described  by  the  classic  thermodynamic  equation 


(Zumdahl,  1989) 


AG  = AG°  + RTln(Q) 


(4-22) 
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where  AG°  is  the  free  energy  change  of  the  system  for  the  given  reaction  with  all  species 
at  1 atmosphere  of  pressure,  AG  is  the  free  energy  change  of  the  system  for  the  given 
reaction  at  the  given  conditions,  R is  the  gas  constant,  and  T is  the  temperature.  Also, 


Q = 


[H  + ]inne,- 
[H  + ]outer 


(4-23) 


is  the  reaction  quotient,  and  [H  + ]inner  and  [H  + ]outer  are  the  hydrogen  ion  concentrations 
of  the  glass-fill  solution  and  the  surrounding  fluid,  respectively.  Substituting  the 
definitions 


AG  :=  -nFE 


(4-24) 


AG°:=-nFE°  (4-25) 

where  n is  the  number  of  electrons  exchanged  in  the  reaction,  F is  the  electrical  charge  of 

a mole  of  electrons,  and  E and  E°  are  the  half-cell  potentials  created  by  the  reactions 

respectively  at  the  current  and  standard  conditions,  gives  the  classical  Nemst  Equation: 


E = E° 


' RT  \ 

JnF  J 


ln(Q) 


(4-26) 


Theoretically,  then,  the  electrical  potential  E-  E°  is  a function  of  T,  n,  and  Q (R 
and  F are  constants).  In  the  glass  electrode  system,  T can  be  considered  a known 
quantity,  as  in  practice  it  is  simply  compensated  for  using  a thermocouple,  and  n is 
always  equal  to  one  (because  hydrogen  ion  always  has  a charge  of  one).  The  remaining 

term  is  the  ratio  of  [H  + ] in  the  glass-fill  solution  ( i.e [H  + ]inner)  to  [H  + ] in  the  solution 


of  unknown  pFl  (i.e.,  [H  + ]outer ).  From  (4-26),  the  natural  log  of  the  ratio  is  proportional 
to  the  electric  potential  being  generated.  A mathematical  transformation  can  be 
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employed  to  change  the  natural  logarithm  to  the  base-ten  logarithm.  Recall  that  pH  is 
defined  as  given  by  (4-1)  as 

pH  = -log[H+] 

Given  that  the  standard  hydrogen  potential  is  the  reference  potential  (i.e.,  E°  :=  0 in  this 
case)  against  which  all  half-reactions  are  assigned  (Zumdahl,  1989),  it  can  be  seen  that 
the  standard  pH  probe  develops  an  electrical  potential  that  is  linear  with  respect  to 

log[H  + ] , i.e.,  with  respect  to  pH.  In  addition  to  reasons  of  brevity,  this  explains  why 
hydrogen  ion  control  is  commonly  referred  to  pH  control. 

4.22.2  Sensing  capabilities 

As  an  inherent  result  of  the  way  that  electrode  sensors  are  designed,  there  are 
numerous  sources  of  error.  Because  of  the  fact  that  two  half-cells  (whether  or  not  they 
are  physically  combined  into  one  probe,  a pH  electrode  always  consists  of  two  separate 
reactions),  a temperature  compensator,  and  an  external  power  supply  are  required,  there 
are  at  least  five  potentials  and  ten  resistances  that  become  involved  in  the  analysis 
(McMillan,  1991).  Furthermore,  although  most  of  these  should  theoretically  be  constant, 
it  is  common  for  them  to  vary  to  some  extent,  resulting  in  a significant  amount  of 
compounded  error.  In  addition,  because  the  sensors  are  constructed  of  glass,  they  are 
generally  limited  to  pressures  of  no  more  than  1 00  pounds  per  square  inch  (psi)  and 

temperatures  of  no  higher  than  100°C,  although  specially  constructed  (and  therefore 
more  expensive)  devices  can  improve  these  bounds  (Moore,  1 990).  In  fact,  by 

reinforcing  the  glass  wall  with  metal,  250  psi  and  140°C  can  safely  be  monitored. 
Unfortunately,  these  sensors  are  extremely  prone  to  erosion  (e.g.  when  exposed  to 
slurries  with  metallic  particles),  corrosion  (hydrofluoric  acid),  and  coating  (black  liquor 
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streams).  Also,  extreme  conditions  (pH  < 3,  pH  > 12,  or  water  content  less  than  40  wt% 
in  the  stream)  tend  to  leach  hydronium  ions  present  in  the  gel  layers,  which  degrade  both 
the  accuracy  and  response  time  of  the  sensors.  When  all  of  these  factors  are  taken  into 

account,  the  best  reliable  accuracy  of  pH  electrodes  that  can  be  expected  is  about  — 0.1 
pH  units,  even  under  well  controlled  laboratory  conditions. 

4.2.2.3  Dynamic  model 

To  formulate  a dynamic  model  for  a pH  electrode,  it  is  easiest  and  most  useful  to 
study  how  one  responds  to  a controlled  change  in  input.  A step  input  can  be  simulated  by 
quickly  transferring  the  pH  probe  from  one  solution  to  another.  An  experimental  curve 
from  such  an  experiment  is  shown  in  Figure  4-3.  In  the  experiment,  the  probe  was  held 


Figure  4-3.  Typical  dynamic  response  of  a pH  electrode  to 
a step  change  in  input. 


115 


in  the  first  solution  of  pH  = 5.7  until  approximately  t = 9 s,  when  it  was  quickly 
transferred  to  a second  solution,  and  the  ion  exchange  was  allowed  to  proceed  to 
equilibrium,  yielding  the  stable  final  pH  reading  of  8.5. 

Based  upon  the  shape  of  the  curve,  a first-order  model  is  appropriate.  A time 
constant  of  approximately  two  seconds  is  observed.  This  is  consistent  with  values 
reported  in  the  literature  (MacMillan,  1 994),  although  it  is  also  stated  that  the  time 
constant  can  vary  (Section  4.3  investigates  the  behavior  of  the  system  when  it  does  vary 
in  the  range  from  1 to  5 seconds).  Figure  4-4  shows  the  model  structure  in  block  diagram 
form 


[Hi 


Figure  4-4.  Block  diagram  of  the  model  of  a traditional  pH 
measuring  electrode. 

where  [H+]  is  the  hydrogen  ion  concentration  in  the  tank,  9elec  is  the  delay  of  the 
electrode,  te]ec  is  the  time  constant  of  the  electrode,  and  pHm  is  the  measured  pH 
produced  by  the  sensor  and  sent  to  the  controller.  Note  that  the  transformation  from  [H+] 
to  pH  occurs  before  the  transfer  function;  this  means  that  a sensor  gain  of  1 
(dimensionless)  is  implied  in  the  transfer  function,  and  is  therefore  omitted  for  clarity. 
The  corresponding  differential  equation  for  the  transfer  function  block  that  represents  the 
probe  model  only  ( i.e the  rightmost  block)  in  Figure  4-4  is  then  given  as 


d(pHm) 


pH, 


PH^-^lec) 


^elec 


Lelec 


dt 


(4-27) 
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It  is  less  certain  what  a reasonable  estimate  of  the  delay  of  the  electrode  would  be,  given 
that  it  is  difficult  to  identify  exactly  at  what  time  the  probe  was  immersed  in  the  second 
solution.  It  is  assumed  that  ion  exchange  begins  quickly  after  the  probe  contacts  the  new 
solution,  creating  little  delay.  A value  of  1 second  is  chosen  throughout  the  study  for  the 
sensor  ( 9elec  = 1 s).  Finally,  since  temperature  could  easily  be  detected  and  passed  into 
the  pH  calculation  algorithm,  it  can  be  assumed  as  a constant  in  simulation.  These 
assumptions  will  be  retained  unless  experimental  results  (see  Section  4.4)  indicate  that 
they  must  be  revised. 

4.2.3  Novel  Optical  Sensors 
4.2.3. 1 Theory  of  operation 

An  alternative  to  pH  sensor  technology  employs  the  use  of  fiber  optics. 
Chromophores  (indicator  dyes  that  are  optically  sensitive,  in  this  case,  to  [H+])  dissolved 
in  a solution  can  be  used  to  accurately  determine  the  pH  of  the  solvent,  with  the  aid  of 
spectrophotometers.  It  has  been  further  demonstrated  that  the  indicator  dyes  can  be 
immobilized  in  a polymer  film  (Seitz,  1988;  Rao  et  al.,  1991),  implying  that  there  is  the 
potential  to  create  sensors  that  can  determine  the  pH  of  flow  streams.  This  is 
accomplished  by  analyzing  the  light  passing  through  a fiber  optic  sensor  coated  with 
polymer  containing  the  indicator  dyes. 

The  apparatus  can  be  placed  into  a flow  stream,  where  the  pH  of  the  fluid 
determines  what  fraction  of  the  immobilized  dye  molecules  accept  a proton.  Because  the 
bond  structures  of  the  protonated  and  unprotonated  molecules  are  different,  their 
absorbance  spectra  also  differ.  An  appropriate  light  source  is  then  placed  near  the  end  of 
the  sensor,  and  the  light  passing  through  is  analyzed  for  absorbency  at  wavelengths 
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critical  to  each  form  of  the  dye.  If  the  total  number  of  sites  is  known,  the  fraction 
protonated,  and  therefore  the  [H+],  can  be  deduced. 

The  dye  or  dyes  that  are  present  in  the  polymer  obey  the  Henderson-Hasselbalch 
equation  (Zumdahl,  1989) 


pH  = pK,-logE^  (4-28) 

[A  ] 

which  establishes  the  relationship  between  pH  and  the  pKa  of  an  acid  as  a function  of  the 
relative  concentration  of  its  protonated  and  unprotonated  forms.  Values  of  pKa  for 
various  acids  are  well  documented,  leaving  the  measurement  of  acid  and  conjugate  base 
concentrations  as  the  primary  obstacle.  Note  that  this  weak  acid  and  its  accompanying 
dissociation  constant  are  not  related  to  acetic  acid  and  its  dissociation  constant.  When 
coupled  with  the  Beer’s  law  (Skoog  et  al.,  1998),  the  spectrophotometer  is  used  to 
accomplish  the  measurement  of  the  needed  concentration.  Beer’s  Law  is 

A = abc  (4-29) 

where  A is  the  absorbance,  a is  the  absorptivity  (a  proportionality  constant),  b is  the  path 

length,  and  c is  the  chromophore  concentration.  By  choosing  an  optimum  wavelength  of 

light,  and  knowing  the  absorptivities  of  the  two  species  that  constitute  the  immobilized 

dye  being  used,  the  thickness  of  the  film  (i.e.,  the  path  length)  and  total  density  of  the 

chromophore  sights  in  the  polymer  (i.e.,  the  concentration),  the  absorbance  detected  can 

be  correlated  to  the  corresponding  acid  and  base  concentrations. 

Using  Beer’s  law  to  calculate  the  hydrogen  ion  concentration  is  a straightforward 
process.  First,  a wavelength  of  light  is  selected  such  that  the  absorptivity  of  cHAopt  is 
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high  and  the  absorptivity  of  cA-  is  low,  or  vice  versa.  Next,  knowing  the  path  length  b, 

opt 

the  absorbance  can  be  expressed 


A = aHAoptbcHAop,+aA  bcA  (4'3°) 

opt  opt 

where  the  only  remaining  unknown  terms  are  A,  cHAo  t , and  cA-  . Since  the  total 

" opt 


chromophore  concentration  is  also  known,  is  observed  that 


cHAopt  =cAop«-cA  (4-3  U 

opt 

in  a fashion  analogous  to,  but  not  directly  related  to,  (4-16).  Finally,  (4-31)  is  substituted 
into  (4-30)  and  the  spectrometer  is  used  to  obtain  a measurement  of  absorbance  A.  After 
suitable  rearrangement 


A aHA0ptbcAopt 

V°Pt  aA  b — aHAoptb 

opt 


(4-32) 


With  cA-  determined,  (4-31)  can  easily  be  used  to  find  cHAo  ( , and  both  can  then  be 

opt  p 


substituted  into  the  final  term  of  (4-28). 

As  with  glass  electrodes,  there  are  physical  limitations  inherent  to  the  design  of 
polymer  sensors  that  must  be  recognized,  analyzed,  and  accounted  for  (Janowiak,  2003). 
First,  the  pKa  of  an  indicator  dye  can  be  affected  when  it  is  immobilized.  In  addition,  it 
has  been  found  that  both  the  ionic  strength  of  the  flow  stream  as  well  as  the  temperature 
have  significant  effects  on  the  apparent  pKa  value  of  the  acid.  Ionic  strength  can  be 
deduced  from  the  conductivity  of  a solution,  which  can  be  measured  easily  by  assembling 
an  electrical  circuit  around  the  solution.  Furthermore,  temperature  is  also  easy  to 
measure,  and  its  dependency  is  no  additional  burden  when  compared  to  traditional  pH 
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sensors.  A model  has  been  proposed  (Chang  et  al.,  1995;  Janowiak,  2003)  to  describe  the 
effects  of  ionic  strength  and  temperature  on  the  observed  pKa  value: 

PKa  = pKa  + P log(Q)  + M ^ (4-33) 

where  pK“  is  simply  a known  reference  value  for  pKa,  P is  an  interaction  factor  between 
the  dye  and  the  electrolytes  present,  Cl  is  the  conductivity  of  the  solution,  and  M is  a 
parameter  related  to  the  activation  energies  of  dye-hydrogen  ion  reactions.  Assuming 
that  this  is  an  accurate  model,  the  pKa  term  in  the  Henderson-Hasselbalch  equation  (4-28) 
can  be  well  defined.  As  in  the  case  of  temperature,  many  sensors  are  available  to  make 
conductivity  measurements.  With  conductivity  and  temperature  measurements,  the  pKa 
of  the  immobilized  dye  can  be  calculated.  Therefore,  using  a spectrophotometer  and 
Beer’s  Law  to  determine  cA  and  cHAo  t , the  ratio  cHAo  t /cA-  can  easily  be  used 

opt  P p opt 

with  (4-28)  and  (4-33)  to  determine  the  pH  of  a solution  in  which  the  probe  is  immersed. 

4.2.3.2  Sensing  capabilities 

As  a result  of  their  simplicity,  optrode  sensors  offer  several  physical  advantages 
to  conventional  probes  (Seitz,  1984;  Janowiak,  2003).  First,  only  one  sensor  is  needed, 
as  opposed  to  two  for  the  conventional  configuration  (recall  that  even  combination  probes 
contain  two  half-cell  reactions);  fewer  signals  limit  the  potential  for  compounded  error. 
Also,  optrodes  are  much  more  rugged  than  their  glass  counterparts,  allowing  higher 
operating  temperatures  and  pressures.  Individually,  they  have  comparably  limited 
dynamic  pH  range  (Seitz,  1984),  but  are  constructed  using  optic  fibers  of  very  small 
diameter  (Borman,  1992),  which  can  potentially  be  clustered  together  to  achieve  equal  or 
greater  dynamic  range  than  conventional  electrodes.  In  addition,  they  are  much  less 
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prone  to  erosion,  corrosion,  and  coating.  Because  the  dyes  can  be  bonded  directly  into 
the  polymer  matrix,  leaching  is  not  a significant  problem,  which  means  that  signal  drift  is 
negligible.  Also,  optrodes  sense  at  the  speed  of  light,  and  are  therefore  limited  by  the 
dynamics  of  the  spectrophotometer  and  dye-hydrogen  ion  reactions,  making  them 
potentially  quicker  than  traditional  electrodes.  Finally,  the  most  attractive  characteristic 

of  the  optical  probe  is  that  pH  measurement  is  linear  with  respect  to  log  — ^ , and 

[A-] 


therefore  the  f H ' 1 measurement  is  linear  with  respect  to  jn  ^is  fashion  the 

[A  ] 

optical  sensor  promises  to  achieve  an  objective  that  may  potentially  outweigh  other 
shortcomings,  namely  linear  feedback.  This  property  could  serve  to  enhance  closed-loop 
control  performance. 

4.2.33  Dynamic  model 

Similar  to  the  electrode,  the  optical  probe  works  via  a concentration-driven  ion 
exchange,  and  it  is  assumed  to  be  described  by  a first  order  model.  In  this  case,  however, 
the  two  species  of  the  immobilized  indicator  dye  are  either  a weak  acid-conjugate  base  or 
weak  base-conjugate  acid  pair.  Similarly  to  the  weak  acid  in  the  influent  stream  (4-5), 
the  total  concentration  and  concentrations  of  protonated  and  unprotonated  species  are 

related  by  a unique  dissociation  constant,  K°pt  as  given  by 


KOpt  , [H+][A-]opt 
3 ' [HA]opt 


(4-34) 
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where  [A  ]opt  is  the  concentration  of  unprotonated  sites  and  [HA]opt  is  the  concentration 
of  protonated  sties.  Note  that  [A~]opt  and  [HA]opt  are  the  same  quantities  given  by  cA- 

opt 

and  cHAopt  in  Section  4.2.3. 1. 

As  with  the  electrode,  the  actual  value  of  the  time  constant  is  assumed  to  vary  for 
different  process  conditions.  Figure  4-5  is  a block  diagram  description  of  the  model 


Figure  4-5.  Block  diagram  of  the  model  of  an  optical  pH 
measuring  probe. 

where  0 is  the  delay  of  the  optrode,  x is  the  time  constant  of  the  optrode, 


[HA] 


opt 


[A"]' 


opt 


IS 


the  real-time  acid-to-base  ratio  of  the  system, 


/r[HA]opt>v 

[A~]opt 


is  the  measured  acid-to-base 


ratio  sent  to  the  controller,  and  [H  + ] is  the  measured  hydrogen  ion  concentration 
measured  by  the  sensor  and  sent  to  the  controller.  Note  that  while  the  sensor  in  fact 
measures  the  acid-base  ratio,  it  is  trivial  to  multiply  this  quantity  by  the  dissociation 

constant  of  the  optrode  K°pt  to  obtain  [H  + ] , and  this  conversion  is  therefore  included  in 
the  model.  As  with  the  electrode  model,  there  is  a sensor  gain  of  1 (dimensionless)  that  is 
omitted  for  clarity.  The  corresponding  differential  equation  for  the  transfer  function 
block  that  represents  the  probe  model  only  (i.e.,  the  middle  block)  in  Figure  4-5  is  then 


given  as 
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[A~]opt  jm  i f[HA]opt 


dt  VUA"]0P' 


In  order  to  compare  the  optrode  performance  directly  to  electrode  performance,  a the 
identical  range  of  time  constant  values  varying  form  1 to  5 seconds  is  selected  for 
simulation.  In  addition,  to  preserve  the  direct  comparison  between  the  two  sensors,  a 
value  of  0opt  = 1 s (recall  that  0e)ec  = 1 s as  well)  is  chosen  throughout  the  simulation 

study.  Finally,  both  temperature  and  conductivity  can  be  assumed  to  be  readily 
measurable  (i.e.,  pK°pt  is  known  exactly)  eliminating  the  need  to  use  (4-33)  in  the 
simulation  model. 


A simulation  study  is  performed  to  examine  if  there  an  advantage  to  using  an 


utilizes  the  plant  model  P,  controller  model  C,  and  measurement  device  (electrode  and 
optrode)  models  M described  in  Section  4.2.  Throughout  the  study,  a number  of 
parameters  remain  unchanged:  namely,  simulation  start  time  is  always  equal  to  t = 0 s, 
step  changes  always  occur  at  t = 20  s,  simulation  end  time  is  always  equal  to  t = 200  s, 
and  each  trial  uses  simple  Euler  algorithm  with  a stepsize  of  0.1  s.  In  addition,  the 
volume  is  defined  as  V = 1 .75  L , and  the  concentration  values  [A]  = [NaOH]  = 0. 1 M 
and  the  flowrate  value  Q2  = 0.0025  L / s are  valid  for  all  initial  steady-state  operating 
points  throughout  the  simulation  study.  Finally,  it  should  be  noted  that  the  equivalence 


4.3  Closed-Loop  Simulations 


optical  probe  rather  than  a traditional  pH  electrode  for  closed-loop  control.  The  study 


point  for  this  system  is  always  pH  = 8.728  in  the  study. 
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Each  closed-loop  trial  is  tuned  separately  to  guarantee  the  optimal  response 
according  to  the  integral  of  time-averaged  error  (ITAE)  criterion  (Seborg,  Edgar,  and 
Mellichamp,  1989;  Ogunnaike  and  Ray,  1994) 


f 

ITAE  = J t e dt  (4-36) 

o 

where  tf  is  the  final  time,  t is  the  time,  and  e is  the  error  between  pH  value  and  pH 
setpoint.  One  major  benefit  of  this  method  is  that  the  best  theoretical  response  attainable 
can  be  identified,  without  additional  complication  introduced  by  the  fact  that  optimal 
tunings  for  controllers  using  either  an  optrode  or  electrode  vary  throughout  the  operating 
range.  In  other  words,  introducing  different  amounts  of  additional  error  by  fixing  the 
controller  tunings  is  avoided.  Finally,  the  ITAE  is  calculated  for  the  period  between  step 
change  ( t = 20  s)  and  the  end  of  the  simulation  ( t = 200  s).  Note  that  a few  simulations 
are  unable  to  reach  final  steady-state  before  t = 200  s;  in  these  cases,  the  simulation 
ending  time  is  increased  for  the  purpose  of  identifying  the  optimal  controller  tunings 
only.  After  the  optimal  tunings  are  identified,  they  are  fixed,  the  simulation  is  re-run  with 
a simulation  stop  time  of  t = 200  s,  and  the  resulting  ITAE  value  recorded. 

A set  of  tuning  heuristics  for  the  definition  of  an  optimal  tuning  of  the  PI 
controller  is  developed.  An  exhaust  numerical  search  is  employed  with  an  adaptive  step 
size  for  gridding,  as  explained  in  Appendix  I.  For  each  of  the  simulation  variations 
described  later  in  this  section,  the  controller  gain  Kc  is  identified  to  two  significant 
figures  (as  defined  in  Appendix  I),  and  the  integral  time  constant  x,  of  the  controller  is 
determined  to  one  significant  figure.  Finally,  while  the  gain  Kc  does  not  need  to  be 
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restricted  in  its  range,  the  restriction  10  s < r,  < 1000  s is  imposed  as  a practical 
limitation  to  unrealistically  extreme  values. 

In  addition,  the  ITAE  criterion  is  applied  to  the  actual  pH  in  all  cases,  rather  than 
the  measured  pH,  which  is  only  used  for  the  feedback  signal.  This  is  particularly 
important  for  the  cases  when  noise  is  injected,  since  additional  error  would  be  artificially 
introduced  by  the  noisy  signal  near  the  end  of  each  simulation.  Finally,  it  should  be 
noted  that  ITAE  applied  to  errors  in  terms  of  the  respective  control  variables  (i.e.,  pH  for 

the  electrode  and  [H  + ] for  the  optrode)  cannot  be  compared  directly.  Therefore,  while 
the  error  for  the  electrode  is  calculated  directly 

eeiec  = PHsp-PH  (4-37) 

where  ee]ec  is  the  error  when  using  an  electrode,  pH  sp  is  the  setpoint  in  units  of  pH  (this 

quantity  is  specified  by  the  user,  and  therefore  always  available),  and  pH  is  the  actual  pH 
of  the  system  (as  opposed  to  the  measured  pH).  In  contrast,  the  error  for  the  optrode  is 
converted  to  pH  units  by  defining  the  error  as 

eoPt=pHsp-(-1ogio([H  + ]))  (4-38) 

where  eopt  is  the  when  using  an  optrode,  and  [H^  ] is  the  actual  hydrogen  ion 

concentration  of  the  system. 

The  control  diagrams  for  both  the  electrode  and  the  optrode  are  derived  from  the 
basic  closed  loop  system  shown  in  Figure  4-1 . The  information  specific  to  the  control 
loop  when  an  electrode  is  used  as  a measurement  device  is  shown  in  Figure  4-6.  Features 
present  in  Figure  4-6  that  are  not  present  in  Figure  4-1  include  disturbances  Q,,  C, , and 
C2  as  well  as  the  injected  noise  signal  n. 
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Figure  4-6.  Block  diagram  of  the  control  loop  when  using 
an  electrode  for  pH  control. 

where  C,  is  equivalent  to  [A]in  and  C2  is  equivalent  to  [Na+]jn  as  given  in  Figure  4-2. 
The  diagram  for  a control  loop  using  an  optrode  as  a measurement  device  is  shown  in 
Figure  4-7,  and  includes  the  additional  features  included  in  Figure  4-6.  Note  that  the 
setpoint  must  pass  through  an  appropriate  function  to  be  in  the  same  units  as  the  feedback 
signal. 

In  each  case,  the  process  begins  with  the  probe  generating  a signal  that  is  a 
function  of  the  hydrogen  ion  concentration  [H  + ] of  the  fluid  in  which  it  is  immersed. 
Noise  n may  then  added  to  the  output  of  the  probe,  and  that  measured  value  is  compared 
to  the  setpoint  by  the  controller.  The  controller  then  determines  the  proper  input  Q2  to 
send  to  the  system.  Finally,  the  system  adopts  a new  hydrogen  ion  concentration  value  as 
a result  of  the  most  current  values  of  Q, , Q2,  C, , and  C2 . 
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Figure  4-7.  Block  diagram  of  the  control  loop  when  using 
an  optrode  for  hydrogen  ion  concentration  control. 

The  total  number  of  simulated  responses  investigated  is  320,  a number  that  results 
from  each  variation  of  six  different  parameters  that  are  allowed  to  be  assigned  more  than 
one  value.  The  six  variables  are:  (/)  sensor  choice,  of  which  there  are  two,  (ii)  value  of 
the  sensor  time  constant,  of  which  there  are  five  (iii)  , whether  or  not  noise  is  injected,  for 
which  there  are  two  choices,  ( iv ) whether  the  response  is  due  to  a servo  or  regulation 
study  (these  terms  are  defined  later  in  this  section),  of  which  there  are  two  choices,  (\;) 
whether  the  disturbance  is  a step  increase  or  a step  decrease,  for  which  there  are  two 
choices,  and  (iv)  one  of  four  initial  operating  points.  These  six  variables  should  be  noted, 
as  they  are  referred  to  frequently  in  later  sections  of  the  study.  Specifically,  the  time 
constant  of  either  sensor  (i.e,  xelec  or  xopt ) can  be  values  of  whole  numbers  from  1 to  5 

seconds,  inclusive.  In  addition,  in  the  case  of  the  servo  responses,  the  initial  pH  values 
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are  6,  7,  8,  and  9,  for  step-up  responses,  and  7,  8,  9,  and  1 0 for  step-down  responses. 

Also,  for  the  case  of  regulation  responses,  the  initial  pH  values  are  7,  8,  9,  and  10,  while 
the  step-increases  and  step-decreases  in  concentration  always  have  a value  of  ten  percent 
of  the  original  value.  Furthermore,  it  is  a primary  objective  that  the  responses  avoid 
saturation  effects  as  much  as  possible,  and  accordingly,  the  chosen  initial  operating  points 
center  around  the  equivalence  point  of  an  acetic  acid-sodium  hydroxide  system 
(pH  =8.728). 

An  automatic  controller  is  typically  required  to  respond  to  two  types  of  changes  in 
the  system  it  is  charged  with  controlling.  One  is  known  as  a servo  response,  and  occurs 
when  the  setpoint  of  the  system  is  changed,  and  the  controller  must  react  properly  to 
move  the  system  fonn  the  original  operating  condition  to  the  new  desired  operating 
condition.  The  second  type  of  response  is  called  a regulation  response,  and  occurs  when 
an  uncontrolled  system  variable  (as  opposed  to  the  input  Q2  of  this  system)  causes  the 
control  variable  (in  this  system,  hydrogen  ion  concentration)  to  change.  Both  types  of 
responses  can  commonly  occur,  necessitating  that  both  types  of  responses  should  be 
studied. 

4.3.1  Noise-Free  Closed-Loop  Performance 

Figure  4-8  depicts  the  time-domain  response  of  a single  combination  of  variables 
(//)-( v/)  as  given  in  Section  4.3  for  both  an  electrode  and  an  optrode.  In  particular,  the 
time  constant  for  each  sensor  (i.e.,  xe!ec  = xopt ) has  a value  of  two  seconds,  and  there  is 

no  noise  injected  into  the  signal.  Additionally,  curves  represent  a servo  response  to  a unit 
step  decrease  beginning  at  pH  = 1 0 . In  the  top  plot,  the  closed-loop  response  of  the 
electrode  is  given  by  a continuous  line,  while  the  closed-loop  response  of  an  optrode  is 
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given  by  a dash-dot  line.  Furthermore,  the  value  of  the  setpoint  is  also  plotted  using  a 
dotted  line.  The  bottom  plot  gives  the  corresponding  input  manipulations  ( i.e .,  the  value 
of  the  reagent  flowrate  Q2 ) implemented  by  the  controller,  and  uses  the  same  line-type 
convention. 


Figure  4-8.  Noise-free  time-domain  servo  responses  of  the 
electrode  (continuous  line)  and  the  optrode  (dash-dot  line) 
sensors  plotted  against  the  setpoint  (dotted  line)  for  a unit 
step  decrease  in  pH  beginning  at  pH  = 1 0 . 

It  is  clear  that  in  this  case,  the  electrode  responds  more  quickly  than  does  the 
optrode.  While  the  two  responses  have  approximately  the  same  settling  time  (i.e.,  the 
time  required  to  reach  each  respective  final  steady-state),  the  rise  time  (i.e.,  the  time  it 
takes  to  reach  the  final  steady-state  for  the  first  time)  is  much  lower  for  the  electrode, 
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causing  additional  error  to  accumulate  for  the  optrode.  In  fact,  this  is  reflected  in  the 
ITAE  values,  which  are  221  for  the  electrode  and  378  for  the  optrode. 

The  closed-loop  performance  defined  by  the  ITAE  method  is  plotted  in  Figure  4-9 
for  forty  of  the  responses  described  in  Section  4.3.  Specifically,  they  are  noise-free  servo 
responses  to  step  increases.  In  other  words,  of  the  six  variables  given  in  Section  4.3,  (/) 
and  (z'z)  are  not  fixed,  (z'z'z)  is  restricted  to  the  case  of  no  noise,  (z'v)  is  restricted  to  the  case 
of  servo  responses,  and  (v)  is  restricted  to  the  case  of  step  increases.  Finally,  each  of  the 
four  subplots  represents  (vz)  one  of  the  four  possible  initial  operating  points.  Within  each 
subplot,  there  are  two  lines  plotted.  The  electrode  is  depicted  by  the  solid  line,  and  the 
optimal  ITAE  value  is  plotted  as  a function  of  the  time  constant  of  the  sensor,  x . Note 
that  for  this  sensor,  x = xe)ec . Similarly,  the  optrode  ITAE  value  is  depicted  by  a dash-dot 
line,  for  which  x = xopt . Note  that  the  line  that  is  relatively  lower  represents  better 
performance. 

It  should  be  noted  that  both  servo  and  regulatory  responses  are  grouped  into  (7) 
steps-up  and  (2)  steps-down  for  clarity.  An  equally  meaningful  arrangement  based  upon 
the  physics  of  the  system  might  be  (a)  changes  that  require  the  controller  to  drive  the 
system  toward  the  equivalence  point,  and  ( b ) changes  that  require  the  controller  to  drive 
the  system  away  for  the  equivalence  point.  In  other  words,  in  general,  some  of  the 
changes  drive  the  neutralization  reaction  toward  a more  sensitive  regime  ( i.e .,  higher 
plant  gain),  while  others  toward  a more  sluggish  regime  (i.e.,  lower  plant  gain).  Note  that 
this  notion  holds  even  for  regulation  responses,  as  the  disturbance  will  serve  to  move  the 
system  either  toward  or  away  form  the  equivalence  point,  and  it  is  the  job  of  the 
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controller  move  it  in  the  opposite  direction  soon  thereafter.  Trends  pertinent  to  this 
phenomenon  are  also  discussed. 

As  shown  in  Figure  4-9,  there  are  operating  regions  where  either  the  electrode  or 
the  optrode  provides  better  control.  Specifically,  in  subplot  (a),  the  optrode  is  more 
effective  and  in  subplot  (c),  the  electrode  is  more  effective.  The  two  measurement 
devices  work  almost  identically  with  respect  to  minimizing  the  ITAE  criterion  in  subplot 
(b).  Finally,  subplot  (d)  represents  a case  where  the  optrode  must  be  significantly  de- 
tuned to  maintain  stability  at  the  original  operating  point  (i.e.,  pH  = 9 ),  thereby  making  it 


x [s]  x [s] 

Figure  4-9.  Noise-free  servo  performance  of  the  electrode 
(continuous  line)  and  the  optrode  (dash-dot  line)  sensors  for 
steps  increases  in  pH  setpoint.  Figure  (a)  shows  a change 
of  pH  from  6 to  7,  (b)  from  7 to  8,  (c)  from  8 to  9,  and  (d) 
from  9 to  10. 
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extremely  sluggish  in  obtaining  the  final  operating  point  (i.e.,  pH  = 10).  Note  that 
subplots  (a)  and  (b)  represent  circumstances  in  which  the  controller  must  move  the 
system  toward  the  equivalence  point.  Finally,  subplots  (c)  and  (d)  are  around  the 
equivalence  point  and  away  from  the  equivalence  point,  respectively. 

The  plots  shown  in  Figure  4-10  represent  a different  set  of  forty  responses  than  do 
the  plots  of  Figure  4-9.  Specifically,  they  encompass  the  cases  where  the  variables  (/)- 
(vz)  of  Section  4.3  are  the  same  as  give  for  Figure  4-9,  except  that  (v)  is  now  restricted  to 
the  case  of  step  decreases. 


Figure  4-10.  Noise-free  servo  performance  of  the  electrode 
(continuous  line)  and  the  optrode  (dash-dot  line)  sensors  for 
steps  decreases  in  pH  setpoint.  Figure  (a)  shows  a change 
of  pH  from  10  to  9,  (b)  from  9 to  8,  (c)  from  8 to  7,  and  (d) 
from  7 to  6. 
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For  the  set  of  step-decrease  servo  responses  given  in  Figure  4- 1 0,  similarly  to  the 
step-increase  trials  reported  in  Figure  4-9,  operating  regions  exist  where  either  probe  is 
more  effective  at  minimizing  the  error.  For  example,  the  electrode  is  significantly  better 
in  the  case  of  subplot  (a),  which  also  corresponds  to  changes  that  move  the  system  toward 
the  equivalence  point.  In  subplots  (b)  and  (c),  however,  there  is  slightly  better 
perfonnance  obtained  when  using  the  optrode,  which  are  cases  of  moving  toward  and 
around  the  equivalence  point,  respectively.  Finally,  they  are  very  comparable  in  subplot 
(d),  where  the  system  is  being  moved  away  from  the  equivalence  point. 

The  eight  comparisons  of  closed-loop  perfonnance  given  collectively  by  Figures 
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Figure  4-11.  Noise- free  regulation  performance  of  the 
electrode  (continuous  line)  and  the  optrode  (dash-dot  line) 
sensors  for  10%  step  increases  in  influent  acid 
concentration.  Figure  (a)  shows  a pH  setpoint  of  10,  (b)  9, 
(c)  8,  and  (d)  7. 


133 


4-9  and  4-10  do  not  seem  to  suggest  that  an  optrode  is  consistently  better  than  an 
electrode  for  servo-response  pH  control.  In  fact,  each  sensor  seems  to  be  more  effective 
under  certain  conditions,  although  there  is  no  obvious  pattern  to  predict  which  one  will 
be. 

Figures  4-1 1 and  4-12  have  plots  analogous  to  the  Figures  4-9  and  4-10,  but  for 
the  case  of  regulation  trials  rather  than  servo  responses  ( i.e .,  variable  (z'v)  is  now  restricted 
to  the  case  of  regulation  responses).  Specifically,  Figure  4-1 1 depicts  the  ITAE  for  ten 
percent  increases  in  influent  acid  concentration  at  pH  setpoints  of  (a)  10,  (b)  9,  (c)  8,  and 
(d)  7.  Furthermore,  Figure  4-12  reports  the  ITAE  for  ten  percent  decreases  in  influent 
acid  concentration  at  pH  setpoints  of  (a)  7,  (b)  8,  (c)  9,  and  (d)  1 0. 

For  the  set  of  step-up  disturbance  rejection  responses  shown  in  Figure  4-11,  the 
results  are  inconclusive,  with  examples  of  better  performance  with  the  electrode,  better 
performance  with  the  optrode,  and  approximately  equal  performance.  For  example, 
subplot  (a)  favors  use  of  the  electrode,  and  is  indicative  of  a change  that  requires  the 
controller  to  move  the  system  away  from  the  equivalence  point.  Next,  subplots  (b)  and 

(c)  give  the  opposite  result,  indicating  that  the  optrode  would  be  a more  effective  choice; 
these  changes  are  around  the  equivalence  point  and  toward  it,  respectively.  Finally,  in 

(d) ,  there  is  no  substantial  difference  in  error  observed  from  using  one  probe  to  the  other, 
with  the  controller  driving  the  system  toward  the  equivalence  point. 

Figure  4-12  shows  the  set  of  steps-down  in  influent  acid  concentration,  and 
similarly  to  the  results  for  step-up  trials,  operating  regions  appear  to  exist  where  the 
electrode  or  optrode  minimizes  the  error  more  effectively,  or  they  perform  approximately 
the  same.  In  subplot  (a),  the  two  sensors  perform  almost  identically  for  a case  where  the 
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Figure  4-12.  Noise- free  regulation  performance  of  the 
electrode  (continuous  line)  and  the  optrode  (dash-dot  line) 
sensors  for  1 0%  step  decreases  in  influent  acid 
concentration.  Figure  (a)  shows  a pH  setpoint  of  7,  (b)  8, 

(c)  9,  and  (d)  10. 

disturbance  moves  the  system  toward  the  equivalence  point.  Next,  subplots  (b)  and  (c) 
clearly  indicate  that  using  the  electrode  is  the  better  choice,  and  represent  changes  toward 
and  around  the  equivalence  point,  respectively.  Finally,  subplot  (d)  indicates  a marginal 
advantage  for  the  optrode,  and  this  for  a case  that  moves  the  system  away  from  the 
equivalence  point. 


Between  Figures  4-1 1 and  4-12,  eight  comparisons  of  closed-loop  performance 


are  given  for  disturbance  rejection,  and  similarly  to  the  servo  tests,  there  is  no  conclusive 
indication  that  an  optrode  performs  better  than  an  electrode  for  pH  control.  Again,  each 
sensor  seems  to  be  more  effective  under  some  conditions,  and  they  are  approximately 
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equal  under  other  conditions.  Also  of  note  is  that  in  none  of  the  sixteen  cases  given  by 
Figures  4-9,  4-10,  4-1 1,  and  4-12  does  one  line  cross  the  other,  representing  a situation 
where  one  sensor  surpasses  the  other  in  performance  as  a function  of  time  constant. 

4.3.2  Closed-Loop  Performance  with  Realistic  Noise 

The  ITAE  trials  performed  in  Section  4.3.2  are  repeated  with  the  addition  of 
measurement  noise  and  reported  in  an  analogous  arrangement.  Therefore,  Figures  4-13 
and  4-14  represent  the  identical  step-increase  and  step-decrease  servo  responses 
previously  described,  and  Figures  4-15  and  4-16  report  the  results  for  step-increase  and 
step-decrease  disturbance  input  concentration.  White  noise  is  selected  for  injection,  with 
magnitude  sufficient  to  ensure  that  three  standard  deviations  are  equal  to  ±0.1  pH  units 
for  the  electrode  and  ± 0.001  pH  units  for  the  optrode,  as  reported  by  (Garcia-Rubio  and 
Clavette,1991).  Finally,  the  electrode  is  depicted  as  a solid  line  in  each  case,  while  the 
optrode  is  a dashed-dot  line. 

The  addition  of  noise  is  equivalent  to  saying  that  the  plots  shown  in  Section  4.3.2 
depict  the  cases  where  variable  (z'z'z)  of  Section  4.3  is  being  restricted  to  the  case  where 
noise  is  being  injected.  Therefore,  coupled  with  the  fact  that  Figures  4-13,  4-14,  4-15, 
and  4-16  otherwise  analogously  represent  the  same  conditions  as  Figures  4-9,  4-10,  4-1 1, 
and  4-12,  the  forty  responses  shown  in  Figure  4-13  are:  (z)  and  (z'z)  are  not  fixed,  (z'z'z')  is 
restricted  to  the  case  that  noise  is  injected,  (z'v)  is  restricted  to  the  case  of  servo  responses, 
and  (v)  is  restricted  to  the  case  of  step  increases.  Finally,  each  of  the  four  subplots 
represents  (vz)  one  of  the  four  possible  initial  operating  points. 
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Figure  4-13.  Servo  performance  of  the  electrode 
(continuous  line)  and  the  optrode  (dash-dot  line)  sensors 
with  realistic  noise  injected  for  steps  increases  in  pH 
setpoint.  Figure  (a)  shows  a change  of  pH  from  6 to  7,  (b) 
from  7 to  8,  (c)  from  8 to  9,  and  (d)  from  9 to  10. 


Figure  4-13  illustrates  the  set  of  results  obtained  for  step-increase  servo  responses. 


Cases  exist  for  operating  regions  where  the  electrode  provides  better  control,  where  the 


optrode  provides  better  control,  and  where  neither  is  significantly  better  than  the  other. 


In  subplot  (a),  for  a case  that  represents  a move  toward  the  equivalence  point,  the  optrode 
is  clearly  superior.  Moving  into  a region  closer  to  the  equivalence  point,  but  still  moving 
toward  it,  subplot  (b)  illustrates  a case  where  the  two  probes  are  comparable.  Next, 
subplot  (c)  centers  around  the  equivalence  point,  and  the  electrode  is  favored.  Finally,  in 
subplot  (d),  the  optrode  suffers  from  instability  at  the  original  operating  point  (/.<?., 


pH  = 9 ) unless  it  is  detuned;  hence  (d)  should  be  ignored. 
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The  forty  optimization  represented  in  Figure  4-14  are  determined  by  using  the 
same  values  of  the  six  variable  ( i)-(vi ) given  in  Section  4.3  used  to  obtain  Figure  4-13, 
except  that  (v)  is  restricted  to  cases  of  step  decreases  only. 

In  subplot  (a)  of  the  set  of  step-decrease  servo  responses  given  in  Figure  4-14, 
which  is  indicative  of  a perturbation  that  drives  the  system  toward  the  equivalence  point, 
the  electrode  is  shown  to  minimize  ITAE  more  effectively.  In  subplots  (b)  and  (c), 
behavior  is  observed  not  previously  seen  in  this  study.  Specifically,  the  time  constant  of 
the  sensors  determine  which  is  the  better  choice.  For  each,  the  optrode  is  indicated  to 
exhibit  better  performance  for  low  values  of  the  time  constant,  but  the  electrode  is 


Figure  4-14  Servo  performance  of  the  electrode 
(continuous  line)  and  the  optrode  (dash-dot  line)  sensors 
with  realistic  noise  injected  for  steps  decreases  in  pH 
setpoint.  Figure  (a)  shows  a change  of  pH  from  1 0 to  9,  (b) 
from  9 to  8,  (c)  from  8 to  7,  and  (d)  from  7 to  6. 
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indicated  at  higher  values  of  the  sensor  time  constant.  Note  that  (b)  is  in  the  region 
around  the  equivalence  point  (pH  = 8.728)  and  (c)  is  a region  where  pH  is  moving  away 
from  the  equivalence  point  for  the  described  simulation.  Finally,  subplot  (d)  suffers  from 
instability  in  much  the  same  way  that  Figure  4-13  does,  rendering  it  inconclusive. 

It  is  observed  that  in  the  presence  of  noise,  certain  conditions  can  result  in 
unstable  behavior  of  the  neutralization  system.  Therefore,  some  trials  must  be  discarded, 
resulting  in  fewer  meaningful  comparisons.  The  six  useful  comparisons  of  closed-loop 
performance  given  collectively  by  Figures  4-13  and  4-14  do  not  suggest  that  an  optrode 
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Figure  4-15.  Regulation  performance  of  the  electrode 
(continuous  line)  and  the  optrode  (dash-dot  line)  sensors 
with  realistic  noise  injected  for  10%  step  increases  in 
influent  acid  concentration.  Figure  (a)  shows  a pH  setpoint 
of  10,  (b)  9,  (c)  8,  and  (d)  7. 
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performs  better  than  an  electrode  for  pH  control.  In  fact,  each  seems  to  be  more  effective 
under  certain  conditions,  although  there  is  no  obvious  pattern  evident  that  can  be  used  to 
predict  which  sensor  should  be  more  effective.  Furthermore,  there  are  two  cases  in  which 
the  sensor  time  constant  plays  a critical  role.  Note  that  in  Figure  4-14  (b)  and  (c),  the 
electrode  is  less  effective  only  in  the  regime  of  fast  sensor  response;  this  can  be  attributed 
to  the  fact  that  measurement  noise  plays  a larger  role  in  performance  degradation  under 
those  conditions.  Figures  4-15  and  4-16  have  plots  analogous  to  the  Figures  4-13  and  4- 
14,  but  for  the  case  of  regulation  trials  rather  than  servo  responses  (i.e.,  variable  (z'v)  is 
now  restricted  to  the  case  of  regulation  responses).  In  other  words,  Figure  4-15  depicts 
the  ITAE  for  ten  percent  step  increases  in  influent  acid  concentration  for  beginning  pH 
values  of  (a)  10,  (b)  9,  (c)  8,  and  (d)  7.  Similarly,  Figure  4-16  reports  the  ITAE  for  ten 
percent  step  decreases  in  influent  acid  concentration  for  beginning  pH  values  of  (a)  7,  (b) 
8,  (c)  9,  and  (d)  10.  Finally,  the  optrode  minimum  error  is  represented  by  the  dotted  line, 
while  the  electrode  minimum  error  is  represented  by  a solid  line. 

The  set  of  step-increase  disturbance  rejection  responses  shown  in  Figure  4-15  represents 
the  first  indication  of  an  obvious  trend  to  this  point  in  the  simulation  study.  First,  subplot 
(a)  must  be  ignored,  since  it  also  suffers  from  instability  problems.  Consequently, 
subplots  (b),  (c),  and  (d)  remain,  and  each  indicates  that  the  optrode  is  better  able  to 
minimize  the  ITAE,  despite  the  fact  that  conditions  moving  both  toward  and  away  from 
the  equivalence  point  are  represented. 

The  set  of  step-decrease  servo  responses  with  noise  are  given  by  Figure  4-16,  and 
contain  examples  that  indicate  that  either  the  optrode  or  electrode  can  perform  better, 
depending  upon  the  conditions.  For  example,  the  neutralization  process  is  moving 
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Figure  4-16.  Regulation  performance  of  the  electrode 
(continuous  line)  and  the  optrode  (dash-dot  line)  sensors 
with  realistic  noise  injected  for  10%  step  decreases  in 
influent  acid  concentration.  Figure  (a)  shows  a pH  setpoint 
of  7,  (b)  8,  (c)  9,  and  (d)  10. 

toward  the  equivalence  point  in  subplot  (a),  and  the  optrode  is  a significantly  better 


choice.  However,  in  subplots  (b)  and  (c),  the  electrode  is  preferable  for  conditions 


representing  a controller  moving  the  system  toward  and  around  the  equivalence  point, 
respectively.  Finally,  (d)  is  ignored  due  to  stability  problems. 


Between  Figures  4-15  and  4-16,  six  meaningful  comparisons  of  closed-loop 
performance  are  given  for  regulation,  and  similarly  to  the  six  meaningful  servo  tests, 
there  is  no  suggestion  that  an  optrode  performs  better  than  an  electrode  for  pH  control, 
even  given  a significantly  lower  level  of  measurement  noise.  Instead,  the  speed  with 
which  the  sensor  responds  (i.e.,  the  sensor  time  constant)  appears  to  be  at  least  as  critical. 
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Finally,  it  is  observed  that  the  injection  of  noise  can  cause  the  optrode  to  exhibit  superior 
performance  to  that  of  the  electrode  in  the  case  of  very  fast  sensor  response  when  it 
would  not  otherwise;  this  behavior  is  a product  of  the  proportionally  larger  effect  noise 
has  on  performance  when  response  characteristics  are  fast. 

4.4  Design  of  an  Experimental  System 
An  important  aspect  of  validating  simulation  results  is  being  able  to  confirm  them 
with  real  data.  To  do  this,  an  experimental  pH  control  is  constructed.  Then  a sensor  is 
wired  to  a computer  with  data  acquisition  capabilities  and  control-capable  software. 
Finally,  the  computer  is  wired  to  the  control  valve  of  the  reagent  stream.  In  addition,  the 
proper  algorithms  are  written  and  implemented  with  the  software. 

There  are  two  important  reasons  for  developing  the  simulation  model  parallel  to 
the  experiment  setup.  First,  no  optical  sensor  is  available.  Therefore,  only  simulations 
are  possible  for  that  portion  of  the  project.  However,  if  experiments  with  an  electrode 
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Figure  4-17.  Flow  diagram  of  a proposed  way  to  compare 
dynamic  response  of  the  sensors  with  no  optical  sensor 
available. 
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match  simulations  using  an  electrode  model,  and  the  models  for  the  sensors  are 
developed  with  equal  care,  extrapolation  can  be  made  with  confidence  about  the  behavior 
of  the  system  under  control  using  an  optical  sensor.  The  process  is  illustrated  as  a 
flowchart  in  Figure  4-17. 

The  second  reason  for  relying  upon  simulations  is  efficiency.  They  are  less 
expensive,  less  dangerous,  and  much  faster.  Troubleshooting  can  be  done  in  simulation, 
greatly  speeding  the  processes  of  controller  tuning  and  studying  system  behavior  at  a 
variety  of  operating  points.  It  is  only  then  that  real  experiments  are  performed  to  confirm 
the  simulation  results. 

4.4.1  Apparatus 

4.4. 1.1  Experimental  hardware 

The  core  of  the  experimental  setup  is  a pre-manufactured  unit  made  by  Armfield 
and  designated  PCT  16  (pH  control  accessory).  A schematic  of  the  equipment  is  shown 
in  Figure  4-18.  The  most  important  components  are  the  sump  tanks,  pumps,  flow  meters, 
control  valve,  impeller,  and  mixing  tank. 

The  process  begins  with  effluent  being  present  in  the  left  tank  and  reagent  in  the 
right  tank.  Pumps  with  recycle  streams  send  a steady  flow  of  each  through  separate  flow 
meters.  The  flow  meter  for  the  effluent  is  adjusted  to  allow  the  desired  flowrate. 
Dissimilarly,  the  flow  meter  for  the  reagent  is  adjusted  to  allow  the  maximum  desired 
flowrate.  The  effluent  then  passes  directly  to  the  mixing  tank.  The  reagent,  however, 
goes  through  the  control  valve  before  it  reaches  the  tank.  The  control  valve  is  capable  of 
receiving  instructions  from  a controller  that  will  further  restrict  the  flow  to  the  desired 
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rate.  Finally,  the  mixture  is  agitated  in  the  tank,  and  overflows  to  a drain  as  more  fluid 
arrives. 


Figure  4-18.  Schematic  of  the  pH  control  experimental 
apparatus. 

4.4.1. 2 Computer  hardware,  data  acquisition,  and  power  supply 

The  computer  used  to  run  the  software  and  send  and  receive  data  is  a Macintosh 
Quadra  700  fitted  with  a PowerPC  upgrade  card.  The  electrical  signals  passing  to  and 
from  the  computer  are  managed  by  a National  Instruments  NB-MIO-16X  data  acquisition 
card  (DAQ).  Finally,  a custom-built  external  power  supply  is  need  for  two  tasks — 
amplifying  the  signal  from  the  DAQ  to  the  control  valve  and  transforming  the  power 
source  for  use  by  the  impeller.  The  physical  setup  is  shown  in  Figure  4-19. 

For  the  chosen  application,  the  output  signals  from  the  controller  (computer)  are 
voltages  sent  to  the  control  valve  of  the  experimental  apparatus.  They  are  generated  by 
the  algorithms  programmed  into  the  software,  and  are  based  upon  the  input  signals.  In 
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Figure  4-19.  Physical  setup  of  experimental  equipment  for 
pH  control. 

this  system,  the  input  signals  to  the  controller  are  readings  sent  from  the  measurement 
device — the  pH  electrode. 

4.4. 1.3  Computer  software 

National  Instruments  also  produces  Labview  (version  3.0),  the  software  used  for 
programming  the  computer  with  the  control  software.  During  construction  and  operation, 
the  program  presents  the  user  with  two  windows  for  each  file.  The  first  step  is  to  select 
the  flow  diagram  screen  (the  window  mostly  used  during  construction)  and  code  in  the 
transfer  function  for  a PI  controller 


u(k)  = Kc  1 


At 

+ — 


(4-39) 


x,  z-l 

where  At  is  the  sampling  interval  in  seconds  (s),  z is  the  complex  z-domain  argument  (a 
shift  operator),  and  Kc  and  T|  are  the  same  as  defined  in  Section  4.2  (Seborg  et  al.,  1989). 
To  put  this  in  a form  immediately  usable  by  Labview  requires  an  input  and  an  output 
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u(k)  = Kc 


(4-40) 
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Applying  the  shift  operator 


] i-l 


(4-41) 


where  k denotes  the  k-th  sample,  e(k)  is  the  error  (input)  signal  at  sample  k,  u(k)  is  the 
signal  calculated  by  the  algorithm  (output)  at  sample  k,  and  other  values  are  defined 
above. 

4.4.2  Measurement  Noise  Analysis 
4.4.2. 1 Signal  assessment 

Measurements  recorded  with  an  electrical  signal  frequently  suffer  from  the  effects 
of  noise.  The  experimental  apparatus  described  in  Section  4.4.1  utilizes  such 
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Figure  4-20.  Noisy  (raw)  pH  signal  originating  from  the 
pH  meter. 
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measurements,  and  a typical  signal  seen  by  the  computer  for  a seemingly  stable  value  of 
pH  displayed  by  the  pH  meter  is  shown  in  Figure  4-20 

It  is,  of  course,  desirable  to  reduce  the  amount  of  noise  in  the  signal.  First,  the 
nature  of  noise  must  be  investigated.  With  the  results,  it  is  possible  to  choose  the  method 
most  appropriate  for  reducing  it.  Because  noise  can  most  effectively  be  filtered  if  it  is 
random,  identification  and  attenuation  of  non-random  signals  is  the  first  step.  Figure  4- 
21  shows  the  power  spectral  density  (found  by  taking  the  Fourier  Transform  of  the  data 
to  convert  it  to  frequency  domain  representation)  of  the  noisy  signal  in  Figure  4-20.  The 
actual  value  of  each  power  is  not  particularly  meaningful;  it  is  the  relative  area  under 
each  spike  that  is  of  importance. 
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Figure  4-21.  Power  spectral  density  of  the  noisy  signal  in 
Figure  4-20. 
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There  are  essentially  two  signals,  one  at  0 hertz  (Hz),  or  cycles  per  second,  and  the  other 
at  approximately  60  hertz.  The  signal  at  0 Hz  is  the  constant  signal  that  is  desired,  while 
the  60  Hz  signal  is  noise  being  introduced.  Because  the  alternating  current  of  the  power 
supply  runs  at  or  near  60  Hz,  this  result  is  quite  plausible.  If  a low-pass  fdter  can  be 
designed  to  eliminate  or  minimize  the  60  Hz  component,  the  random  nature  of  the 
remaining  signal  can  be  confidently  assumed.  There  are  two  main  categories  of  filters — 
hardware  and  software,  and  each  are  explained  further  in  turn. 

4.4.2.2  Hardware  filtering 

Before  the  signal  ever  gets  to  the  computer,  it  is  possible  to  condition  the  signal 
using  basic  electronic  components — resistors,  capacitors,  etc.  While  more  complicated 
schemes  are  possible,  a traditional  RC  (resistor-capacitor)  circuit  is  chosen  for  its 
simplicity  and  ease  of  analysis.  Essentially,  they  are  low-pass  filters  with  a cutoff 
frequency  equal  to  RC.  The  unit  of  R is  ohms,  and  the  unit  of  C is  Farads,  the  product  of 
which  is  frequency  in  hertz.  The  physical  configuration  of  an  RC  circuit  is  illustrated  in 
Figure  4-22.  The  attenuation  rate  of  the  filter  is  smooth  and  continuous,  so  the  60  Hz 
signal  is  not  completely  eliminated.  However,  it  drops  off  sharply  at  higher  frequencies, 
so  appropriate  selection  of  the  term  RC  can  ensure  adequate  performance.  Additionally, 
because  there  is  no  voltage  (signal)  drop  across  the  filter,  multiple  filters  in  series  can  be 
employed  without  adverse  effects.  In  this  case,  a cutoff  frequency  of  10  hertz  is  selected. 


R 


Figure  4-22.  Schematic  representation  of  an  RC  circuit. 
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and  two  such  filters  are  used. 

Although  the  signal  appears  unchanged  to  the  naked  eye  in  the  time  domain, 
Figure  4-23  proves  that  the  contribution  of  the  60  Hz  signal  is  greatly  reduced. 

Using  an  RC  circuit  has  both  advantages  and  disadvantages.  First,  they  are  of 
negligible  cost,  and  extremely  easy  to  implement.  Second,  they  do  not  introduce  any 
time  delay  to  the  system.  However,  they  must  appear  before  the  DAQ  component, 
meaning  there  will  always  be  the  potential  for  some  signal  corruption  after  it  has  been 
filtered.  Finally,  the  discharging  of  the  capacitor  causes  spikes  in  the  reading  (refer  to 
Figure  4-3). 


Figure  4-23.  Power  spectral  density  of  the  noisy  signal  in 
Figure  4-20  after  RC  filtering. 
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4.4.2.3  Software  filtering 

After  the  signal  is  delivered  to  the  computer,  there  are  additional  ways  to 
condition  the  signal  that  is  ultimately  sent  to  the  controller.  They  are  software  filters,  and 
two  of  the  most  basic  types  are  considered — exponential  and  moving  average.  They  are 
similar  in  concept,  using  some  type  of  arithmetic  average  of  a predetermined  number  of 
data  points.  When  the  associated  parameters  are  chosen  carefully,  the  result  is  a 
significantly  smoothed  version  of  the  incoming  data  stream.  See  Appendix  J for  the 
equations  of  each,  as  well  as  qualitative  visual  representations  of  how  they  affect  a raw 
set  of  data. 

Using  algorithms  in  the  control  software  can  help  alleviate  some  of  the 
shortcomings  of  hardware  filtering.  Because  the  signal  has  already  been  digitized, 
additional  noise  further  downstream  is  not  a consideration.  In  addition  (for  the  same 
reason),  the  spikes  in  a signal  that  has  been  filtered  continuously  can  be  smoothed. 
Unfortunately,  even  though  they  are  also  inexpensive  and  easy  to  use,  they  introduce  the 
most  critical  problem  not  present  with  a hardware  filter — delay.  There  is  no  way  to 
circumvent  it,  since  each  must  spend  a finite  amount  of  time  waiting  on  sufficient  data  to 
arrive  before  it  can  make  the  calculation.  As  a result,  selection  of  the  filter  parameter 
(each  has  one)  represents  a tradeoff  between  smoothness  and  introduced  delay. 

The  final  decision  of  which  type  to  use  is  made  from  experience,  when 
experiments  on  the  actual  system  indicate  which  type  tends  to  achieve  the  best 
smoothness/delay  compromise.  In  this  case,  preliminary  results  appear  to  show  that  a 
moving  average  filter  will  work  more  effectively;  however,  this  is  open  to  scrutiny,  and 
could  be  amended  as  the  experimental  design  evolves. 
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4.5  Conclusions 

The  simulation  portion  of  the  work  is  complete.  A representative  group  of  both 
servo  responses  and  disturbance  rejection  responses  have  been  performed  for  both  the 
novel  optical  pH  sensor  and  the  traditional  pH  electrode.  Furthermore,  each  was 
simulated  for  multiple  values  of  sensor  time  constant.  Tuning  was  done  for  each 
individual  case  via  trial-and-error  to  ensure  the  best  settings  for  that  particular  trial,  based 
upon  ITAE  as  the  performance  criterion. 

The  optical  sensor  is  not  consistently  observed  to  guarantee  better  closed-loop 
performance  when  compared  to  the  traditional  pH  electrode.  In  fact,  there  is  a significant 
number  of  cases  where  the  electrode  is  a better  choice.  While  it  is  generally  accepted  that 
an  entirely  linear  system  (i.e.,  controller,  plant  and  sensor)  is  the  easiest  to  control,  the 
plant  in  this  case  is  highly  non-linear.  The  results  to  this  point  support  the  speculative 
idea  that  a linear  sensor  will  not  always  provide  the  best  control  performance  for  a non- 
linear plant.  Rather,  a non-linear  sensor  that  is  non-linear  in  precisely  the  proper  way 
may  serve  to  somewhat  ‘cancel’  the  non-linearity  of  the  plant  under  certain 
circumstances,  thus  creating  a more  linear  overall  system.  This  hypothesis,  if  proven 
correct,  could  theoretically  be  applied  to  other,  unrelated,  non-linear  systems  in  order  to 
improve  performance. 

Another  possibility  that  may  explain  the  results  is  that  pH  units  are  the  basis  of 
reference  for  the  error  criteria,  as  discussed  for  (4-38).  These  units  were  chosen  because 
most  rules  concerning  the  allowable  hydrogen  ion  concentration  in  a process  stream  are 
written  in  terms  of  pH;  this  is  natural  given  that  pH  is  the  units  being  measured  in  such 
streams.  However,  the  true  issue  is  whether  one  hundred  times  the  desired  hydrogen  ion 
concentration  is  really  only  twice  as  large  an  error  as  ten  times  the  desired  hydrogen  ion 
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concentration  (i.e.,  in  ITAE,  an  error  of  two  pH  units  is  weighed  twice  what  an  error  of 
one  pH  unit  is  weighed  for  the  same  instant  in  time).  If  quality  standards  were  revised  to 
be  expressed  in  the  true  unit  of  interest — hydrogen  ion  concentration — a similar  study 
might  reveal  a clearer  indication  that  the  fiber-optic  linear  sensor  is  indeed  consistently  a 
better  choice  for  closed-loop  pH  control  than  is  a pH  electrode. 

Despite  the  lack  of  obvious  superiority  for  either  sensor,  a number  of  conclusions 
can  be  made  about  the  linear,  optical  sensor  and  its  potential.  First,  the  time  constant 
plays  a critical  role  with  respect  to  performance  for  both  sensors,  so  the  optical  sensor  can 
yield  a clear  advantage  to  the  electrode  if  it  can  be  designed  to  respond  quickly.  Also,  at 
the  very  least,  an  optrode  whose  output  is  converted  to  pH  units  would  perform  better 
than  an  electrode  simply  due  to  the  lower  measurement  noise  associated  with  it.  Finally, 
the  tunings  for  the  electrode  vary  far  more  over  a given  operating  range  than  do  those  of 
the  optical  probe;  for  applications  where  the  controller  tunings  cannot  be  adapted  in  real- 
time efficiently,  the  fiber-optic  sensor  will  provide  better  control  as  the  operating  point 
changes  dynamically. 


CHAPTER  5 

CONCLUSIONS  AND  FUTURE  WORK 
A wastewater  disinfection  reactor  that  utilizes  chlorine  has  been  realistically 
modeled.  A standard  cascade  control  system  is  recommended  to  reject  the  effect  that 
ammonia  demand  has  on  residual  chlorine  levels.  However,  many  possible  master 
controller  algorithms  fail  to  deliver  the  desired  outlet  chlorine  concentration  when  the 
flow  rate  through  the  plant  varies  significantly.  It  is  proposed  that  a supplemental  control 
loop  is  required,  where  a residence  time  controller  (RTC)  manipulates  the  volume  of  the 
reactor,  thereby  dramatically  increasing  the  effectiveness  of  the  cascade  system.  Finally, 
the  RTC  can  be  implemented  in  either  feedforward  or  feedback  form,  and  the  appropriate 
choice  is  application-specific. 

Better  control  behavior  and  a minimum  average  25-30%  annual  reduction  in  the 
amount  of  chlorine  used  can  be  expected  by  implementing  this  technology.  For  an  8 
million  gallon-per-day  (8MGD)  plant,  the  dynamic  weir  system  can  save  the  treatment 
plant  approximately  32  US  tons  per  year  over  the  fixed-weir  design.  And  this  is  the 
potential  from  flow  equalization  and  a conservative  residence  time  set  point  (3600  s) 
alone.  Florida  regulations,  for  example,  require  only  900  s of  contact  time,  promising  the 
possibility  of  more  chlorine  savings  with  a lower  RTC  set  point.  Finally,  since  much  less 
outlet  chlorine  concentration  fluctuation  is  seen  with  the  dynamic  weir  implementation, 
the  set  point  (i.e.,  presently  2 ppm)  can  be  moved  closer  to  the  minimum  stipulated  by  the 
regulation  (i.e.,  1 ppm  in  this  case).  Taken  together,  the  aforementioned  strategies  can 
accomplish  an  even  more  dramatic  chlorine  savings.  For  example,  for  the  same  8 MGD 
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plant,  if  the  outlet  concentration  set  point  is  fixed  at  1 .2  ppm  and  the  residence  time  set 
point  is  fixed  at  1 800  s,  the  yearly  total  chlorine  would  then  becomes  approximately  52 
US  tons,  representing  a savings  of  up  to  45%.  This  is  possible  because  the  volume  of  the 
chlorine  contact  basin  is  adjusted  in  real-time,  and  the  chlorine  dosage  system  is  always 
operated  optimally.  In  other  words,  wastewater  is  generally  not  retained  within  the 
reactor  longer  than  necessary,  consuming  excessive  chlorine,  as  is  the  case  with  a fixed- 
weir  structure.  Furthermore,  high  chlorine  dosages  that  result  in  the  formation  of 
carcinogens  can  also  be  eliminated.  Finally,  a plant  that  uses  less  chlorine  has  a 
proportionally  lower  chlorine  storage  requirement,  which  helps  promote  a safer 
environment  at  the  treatment  facility. 

Future  work  in  this  area  of  research  would  certainly  involve  experimental  trials. 
To  do  so,  the  permission  of  regulatory  bodies  such  as  the  EPA  would  be  the  first  step. 
Next,  a suitable  wastewater  treatment  plant  would  need  to  be  identified,  and  the  moving 
weir  installed.  The  effectiveness  of  the  new  control  strategy  could  then  be  evaluated.  It 
should  be  noted  that  pilot-scale  trials  are  not  recommended,  since  the  scalability  of  a key 
dimensionless  parameter  of  the  system,  the  Froude  number,  would  be  difficult  to  preserve 
through  any  significant  reduction  in  reactor  size. 

The  derivation  of  run-to-run  proportional-integral  controllers  is  presented.  With 
knowledge  that  robustness  is  described  by  a gain  margin,  a set  of  tuning  heuristics  is 
given  for  static  plants  with  significant  delay.  With  them,  stable  plant  behavior  can  be 
guaranteed.  In  contrast,  when  a Smith  Predictor  is  used  to  attempt  to  improve 
performance  for  this  type  of  system,  instability  is  observed  even  for  small  modeling  error. 
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When  a Smith  Predictor  is  used,  it  is  observed  that  this  type  of  system  can 
actually  go  from  an  unstable  state  to  a stable  one  (and  back  again)  as  the  difference 
between  true  measurement  error  and  the  estimated  error  of  the  Smith  Predictor  increases. 
Therefore,  it  is  concluded  that  at  least  one  root  of  the  characteristic  equation  of  the 
system  with  a Smith  Predictor  migrates  in  and  out  of  the  stable  region  for  monotonically 
increasing  error  in  delay.  More  importantly,  even  a mismatch  in  delay  of  only  one  unit 
can  cause  instability.  Therefore,  unless  the  delay  of  the  system  is  known  with  certainty  in 
real  time,  it  is  recommended  that  simple  PI  control  be  used  for  run-to-run  static  systems 
with  significant  delay. 

To  further  develop  this  work,  a comprehensive  study  of  the  behavior  of  the  roots 
of  the  characteristic  polynomial  may  be  helpful  in  finding  ways  to  modify  the  controller 
to  be  more  robust  with  a Smith  Predictor.  Without  that  development,  it  is  necessary  to 
use  more  traditional  control  approaches,  even  though  performance  will  be  limited.  A 
possible  starting  point  may  be  to  study  the  discrete  analog  of  the  findings  of  (Adam  et  al., 
2000),  who  have  shown  that  discrete  stability  ranges  exist  for  a Smith  Predictor  applied 
to  a continuous  system  with  an  uncertain  time-delay  parameter.  In  their  work,  the 
frequency  response  in  the  Nyquist  plane  is  shown  to  rotate  about  the  origin  as  the  delay 
value  varies,  which  can  alter  the  number  of  encirclements  of  the  critical  point  (Latchman 
and  Crisalle,  1995),  which  is  a parameter  that  is  central  to  the  stability  of  a continuous 
system. 

The  closed-loop  dynamic  behavior  of  a CSTR  for  effluent  pH  control  has  been 
modeled  with  both  a traditional  pH  electrode  sensor  and  a fiber-optic  sensor  with  a linear 
response  to  hydrogen  ion  concentration.  Furthermore,  a representative  group  of  both 
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possible  responses  were  performed  over  a realistic  range  of  sensor  time  constant.  Each 
trial  was  tuned  individually  based  upon  the  ITAE  performance  criterion.  When 
compared  for  the  same  sets  of  operating  conditions,  there  is  no  indication  that  the  linear 
optical  sensor  performs  better  than  the  traditional  pH  electrode.  In  addition,  there  are 
certain  circumstances  where  the  electrode  is  clearly  more  effective. 

It  is  widely  accepted  that  the  neutralization  process  is  highly  non-linear,  and  the 
results  suggest  that  substituting  a linear  measurement  device  into  the  control  loop  is,  at 
best,  only  a partial  remedy  for  the  control-degrading  non-linear  nature  of  the  system.  In 
fact,  a non-linear  sensor  with  just  the  proper  characteristics  may  better  produce  a quasi- 
linear  composite  system,  although  describing  the  necessary  characteristics  precisely 
would  be  challenging. 

An  additional  explanation  for  the  result  that  the  optical  sensor  is  not  clearly 
superior  to  an  electrode  may  come  from  the  fact  that  their  performance  is  judged  in  terms 
of  pH  units.  If  error  were  calculated  in  terms  of  hydrogen  ion  concentration,  and  the 
loops  retuned  accordingly,  more  obvious  performance  trends  may  manifest  themselves. 
Regardless,  there  are  a few  clear  conclusions  that  can  be  made  about  the  behavior  of  the 
fiber-optic  sensor  with  respect  to  a traditional  electrode.  For  example,  the  time  constant 
is  critical  to  one  measurement  device  outperforming  the  other,  so  a sufficiently  fast 
response  characteristic  would  be  required  to  choose  one  over  the  other.  Also,  it  may  be 
possible  to  simply  convert  the  output  of  an  optrode  to  pH  units,  thus  exploiting  the 
advantage  of  lower  measurement  noise,  with  no  other  source  of  difference  in  their 


respective  responses. 
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Future  work  for  this  area  of  investigation  may  involve  changing  the  definition  of 
the  error  for  the  purpose  of  ITAE  calculation.  In  terms  of  hydrogen  ion  concentration, 
one  sensor  or  the  other  may  prove  itself  to  be  clearly  superior  in  closed-loop  performance 
tests.  Additionally,  experimental  trials  should  be  performed  to  verify  that  the  closed-loop 
system  behaves  as  predicted  for  each  of  the  sensors.  With  experimental  results,  the 
conclusions  drawn  from  the  simulation  study  would  be  validated,  and  additional  ways  to 
exploit  the  nature  of  linear  fiber-optic  pH  sensors  could  be  investigated. 


APPENDIX  A 

DISTRIBUTED  PARAMETER  EQUATION  DERIVATION 
This  appendix  presents  a detailed  derivation  of  the  PDE  that  describes  a PFR 
given  in  Section  2.2.4.  It  begins  with  a conservation  of  mass  balance,  and  uses  the  result 
in  a shell  balance  around  a small  portion  of  the  reactor.  Finally,  it  takes  the  limit  as  the 
incremental  volume  approaches  zero.  By  the  very  derivation  that  results  in  a PDE,  the 
reverse  process  (use  of  difference  equations  to  approximate  one  differential  term)  is 
concluded  to  be  a plausible  alternative  expression. 

Consider  a small  section  of  the  reactor 


Mass  Balance 

Beginning  with  an  overall  mass  balance 
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► 

► 
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d(pAVca) 

dt 


pQz  PQz+Az 


Using  the  fact  that  input  flow  rate  is  equal  to  output  flow  rate  for  this  system 


(A-l) 


Qz  Qz+Az 
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Substituting  (A-2)  into  (A-l)  and  simplifying 


d(pAVca) 

dt 


pAVca  = constant 

Noting  that  water  is  incompressible 


Shell  Balance 


p = constant 
AVca  = constant 


Performing  a shell  balance  around  a small  portion  of  the  reactor 


d[CAVca] 

dt 
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AVca  = AAz 

where  A is  the  cross-sectional  area 


— AAz  = Q[CZ  - Cz+4z]  - k„„CAAVca 
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where  v is  the  velocity  of  the  fluid.  Taking  Hm  yields  a PDE 
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Because  the  equation  is  based  upon  the  differential  volume  model,  it  is  logical  that  it  can 


be  approximated  by  a large  number  of  extremely  small  reactors. 


APPENDIX  B 

PERFORMANCE  ASSESSMENT  OF  THE  STANDARD 
FEEDBACK  CONTROL  CONFIGURATION 

This  appendix  contains  information  related  to  simulations  using  standard 
feedback  control,  including  the  block  diagram,  open-loop  response  of  the  plant,  controller 
tunings,  and  closed-loop  response.  It  is  concluded  in  the  text  of  the  report  that  cascade 
control  should  be  more  effective,  and  the  information  below  bears  this  out.  As  a result, 
the  details  related  to  standard  feedback  are  omitted  from  the  body  of  the  text,  and  are 
presented  below. 

The  block  diagram  for  a standard  feedback  system  (Figure  B-l)  is  very  similar  to 
the  one  given  for  a cascade  feedback  system  shown  in  Figure  2-6  of  the  main  text.  All 
variables  given  are  the  same  as  those  listed  for  the  cascade  system  in  Section  2.4 


w 

Q(t) 
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Figure  B-l.  Block  diagram  of  a wastewater  treatment  plant 
under  standard  feedback  control  (in  addition  to  a ratio 
controller). 

When  the  derived  models  of  the  dosing  area  and  contact  are  used  for  Gda  and  Gca, 


and  the  controllers  are  turned  off,  the  following  open-loop  response  is  obtained 
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Figure  B-2.  Open-loop  response  of  the  entire  plant  when 

subjected  to  a step  change  in  dosing  rate  from  2.926 
3 3 

m /s  to  3.326  m /s  at  time  2 hr. 

The  response  shown  in  Figure  B-2  is  then  fit  to  a FODT  model  exactly  as  those 


in  Sections  2.4.1  and  2.4.2,  and  the  resulting  tuning  parameters  (also  found  through  trial- 
and-error)  are  listed  in  Table  B-l 


Table  B-l.  Estimated  transfer  function  parameters  for 
the  entire  plant  and  resulting  master  controller 
parameters  found  by  trial-and-error  tuning. 


Transfer  Function 

Numerical  Value 

Controller 

Numerical  Value 

| Parameters 

Parameters 

* 

Kca 

1.16  ppm-s/g 

Kmc 

0.0003  g/ppm-s 

j * 

s Tca 

850  s 

Tmc 

2000  s 

0ca 

6430  s 

It  happens  that  these  are  the  tunings  found  for  both  a PI  master  controller  and  a PI 


master  controller  with  a Smith  Predictor.  This  in  unusual,  since  a Smith  Predictor  usually 
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allows  the  operator  to  make  the  tuning  parameters  more  aggressive;  verification  of  this 
result  is  addressed  in  Section  2.6.  Using  these  tuning  parameters,  the  closed-loop 
performance  of  the  standard  feedback  scheme  is  compared  to  that  of  a cascade  scheme 
(using  parameters  given  in  Section  2.4.2)  for  a sinusoidal  disturbance  in  ammonia 

chlorine  demand,  5(t).  The  results  are  given  in  Figure  B-3;  note  that  Q(t)  was  held 

3 

constant  at  0.4  m /s. 


Figure  B-3.  Closed-loop  comparison  of  a feedback  control 
scheme  versus  a cascade  control  scheme  for  a sinusoidal 
disturbance  in  ammonia  chlorine  demand,  5(t). 

Furthermore,  the  response  of  the  system  under  feedback  control  is  virtually  the 
same  regardless  of  what  type  of  master  controller  is  used.  It  is  clear  in  Figure  B-3  that, 
while  the  feedback  control  is  not  bad,  the  cascade  control  is  clearly  superior.  It  is  for  this 
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reason  that  it  has  been  eliminated  from  consideration  in  the  study  as  a viable  control 
scheme. 


APPENDIX  C 

REVIEW  OF  THE  OPERATION  OF  A DISCRETE  CONTROLLER 
APPLIED  TO  A CONTINUOUS  SYSTEM 

Figure  C-l  illustrates  how  a continuous  output  signal  in  the  first  plot  is 


sampled  every  At  time  units.  That  value  is  passed  to  the  digital  control  algorithm  (3-6), 


which  requires  s (note  that  this  s is  different  from  a(k),  the  error  of  batch  k)  time  units  to 


calculate  and  implement  the  new  input.  A zero-order  hold  is  then  used  to  continue  the 


input  until  a new  value  is  calculated  and  implemented  at  the  end  of  the  next  iteration. 


i i 


Figure  C-l.  Graphical  representation  of  the 
implementation  delay  of  a discrete  controller. 

Clearly,  it  is  important  for  e « At  to  ensure  that  the  approximation:  u(k  + e)  is 
valid;  given  the  speed  of  a typical  digital  computer  this  condition  is  easy  to  guarantee. 
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APPENDIX  D 

CONDITIONS  FOR  WHICH  THE  TRANSFER  FUNCTION  (3-8)  HAS  UNSTABLE 

POLE-ZERO  CANCELLATIONS 

The  purpose  of  this  appendix  is  to  show  that  there  are  no  unstable  pole-zero 
cancellations  for  the  transfer  function 


K, 


1 + Kr 


A 

U 


+ R, 


z- 1 


KP 


= 0 


+ R, 


Kpz 


-D 


vz  z — 1 

which  is  shown  as  (3-8)  in  the  main  body  of  the  text.  Multiplying  through  by  the  term 


z(z  - 1 l>  and  rearranging  terms  yields  the  rational  function 


KCKP(R,  + l)zD+1  - KcKpzD 


zD+2  -zD+1  +KCKP(R,  +l)z-KcKp 


= 0 


(D-l) 


Further  factoring  out  the  term  KCKP(R,+1)  gives 


KCKP(R,  +1)- 


,z  R,+r 


zD+2  - zD+1  + KCKP(R,  + l)z-  KCKP 


= 0 


(D-2) 


Recognizing  that  R,  > 0 , it  follows  that  there  may  be  a zero  on  the  unit  circle  only  when 


R,  = 0. 


Consider  the  cases  (i)  R,  > 0 , and  (ii)  R,  = 0 


(i)  If  R,  > 0,  then  0 < — - — <1  and  the  zero  of  (D-l)  is  strictly  inside  the  unit 

R.  + 1 

circle.  This  zero  cannot  therefore  be  cancelled  with  an  unstable  pole. 
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(ii)  If  R,  = 0 , there  is  a zero  of  (D-l)  at  z = 1 , exactly  on  the  unit  circle.  Now  let 
P(z)  = zD+2  -zD+l  + KcKp(R,  +l)z-KcKp  and  after  recognizing  that 
P(l)  = 0 when  R,  = 0 , it  follows  z = 1 is  also  a pole  of  (D-l).  Therefore,  it 
can  be  concluded  that  there  is  always  an  unstable  pole-zero  cancellation  if 
R,  = 0 . Checking  for  repeated  poles  at  z = 1 , divide  P(z)  by  (z  = 1)  to  obtain 


P'(z) 


zD+2-zD+1  . KCKP(R,+1)  KcKp 


z- 1 


+ - 


z — 1 


z- 


z- 1 


(D-3) 


Evaluating  at  R,  = 0 yields 


D+2  _ D+l  k K 

P'(z)  = + Kr  K„  c p 


z- 1 


e p 


z — 1 z — 1 


(D-4) 


or 

P'(z)  = zD+1+KcKp  (D-5) 

Now  since  P'(l)  = 1 + KcKp , it  follows  that  z = 1 is  a duplicated  zero  P(z) 
when  R,  = 0 if  and  only  if  KcKp  = -1 . However,  the  first  necessary 
condition  of  the  Jury  Test  strictly  requires  that 

|kckp|  < i 

so  the  Jury  Test  will  always  determine  whether  or  not  the  system  under 
consideration  is  stable  provide  that  R,  > 0 . 


APPENDIX  E 

VALUES  OF  a AND  (3  FOR  THE  2-PARAMETER  CURVE  FITS  OF  THE 
FUNCTIONS  GENERATED  BY  THE  JURY  TEST 


Delay 

Alpha 

Beta 

0 

0.0700 

89.2029 

1 

0.0800 

85.3576 

2 

0.1100 

67.1585 

3 

0.2239 

25.1400 

4 

0.3502 

14.874 

5 

0.7502 

4.6984 

6 

0.9003 

4.2589 

7 

1.106 

3.6091 

8 

1.3046 

3.2504 

9 

1.7642 

2.2637 

10 

2.0646 

2.0139 

11 

2.2345 

2.0286 

12 

2.4063 

2.0380 

13 

2.5717 

2.0566 

14 

2.7430 

2.0617 

15 

2.9119 

2.0703 

16 

3.0820 

2.0784 

17 

3.2550 

2.0800 

18 

3.4226 

2.0880 

19 

3.6014 

2.0831 

20 

3.7728 

2.0875 
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APPENDIX  F 

CONDITIONS  FOR  WHICH  THE  TRANSFER  FUNCTION  (3-38)  HAS  UNSTABLE 

POLE-ZERO  CANCELLATIONS 

The  purpose  of  this  appendix  is  to  show  that  there  are  no  unstable  pole-zero 
cancellations  for  the  transfer  function  (3-38)  of  a SISO  static  run-to-run  system  with 
measurement  delay  and  a Smith  Predictor.  Equation  (3-38)  shows  that  the  numerator  of 
the  transfer  function  is  the  same  as  for  the  case  when  no  Smith  Predictor  is  employed. 
Therefore,  according  to  Appendix  D,  we  must  find  conditions  under  which  z = 1 is  a root 
of  the  characteristic  equation  when  R,  = 0 . Beginning  with  the  characteristic  equation 
(3-40) 


zD"+D+2  +[KcKp(R,  + 1)  - l]zD+D+1  -(KcK°p)zD"+D-[KcK°(R,  +l)]zD+I 
+(KcKp )zD  + [KCKP (R , + 1 )]zD"+l  - (KCKP )zD"  = 0 
set  R,  = 0 , factor  out  (z- 1) , and  then  set  z = 1 in  all  remaining  terms 


(l  + KcKp)-  Kt.Kp  + KCKP  = 0 (F-l) 

Solving  this  expression  for  zeros  yields 

1 + KCKP  = 0 (F-2) 

Applying  the  first  condition  (3-17)  of  the  Jury  Test  to  (3-40)  shows  that  the  root 

KCKP  = -1  is  forbidden.  There  are  there  possible  cases  of  the  value  of  an , depending 
upon  the  value  of  D.  Consider  the  cases  (i)  D = D°,  (ii)  D < D° , and  (iii)  D>D°. 

(i)  If  D = D° , then  the  final  four  terms  cancel  each  other,  and  the  third  term  in 
(3-40)  becomes  that  with  the  lowest  power  of  z,  therefore  its  coefficient 
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( -KCKP  ) is  by  definition  an . 

(ii)  If  D < D° , then  the  fifth  term  in  (3-40)  becomes  that  with  the  lowest  power  of 
z,  therefore  its  coefficient  ( KCKP ) is  by  definition  an . 

(iii)  If  D > D° , then  the  seventh  term  in  (3-40)  becomes  that  with  the  lowest  power 
of  z,  therefore  its  coefficient  (-KCKP)  is  by  definition  an . 

For  any  of  the  cases,  the  necessary  condition  (3-17)  states  that  |KcKp|  < 1 . Allowing  the 

statement  (F-2)  to  be  true  would  violate  the  condition,  proving  that  there  are  no  unstable 
pole-zero  cancellations  that  could  alter  the  conclusions  of  testing  the  roots  of  the  char- 
cteristic  equation  for  location  relative  to  the  unit  circle. 


APPENDIX  G 

REVIEW  OF  WEAK  ACID-STRONG  BASE  CHEMISTRY 
By  definition,  when  a weak  acid  or  base  is  dissolved  in  water,  its  molecules 
always  exist  in  (at  least)  two  states,  protonated  and  unprotonated.  (Molecules  that  can 
donate  or  accept  multiple  protons  have  additional,  partially  protonated  states.)  The 
relative  abundance  of  the  two  states  is  related  through  the  dissociation  constant,  K.  For  a 
weak  acid,  it  is  designated  Ka  (for  a weak  base,  Kb),  and  is  defined  the  equation 


[H+][A-] 

[HA] 


(G-l) 


where  Ka  is  a constant  related  to  a particular  acid,  [HA]  is  that  acid,  [A  ] is  the  conjugate 
base  of  that  acid,  and  [H  ] are  free  existing  hydrogen  ions;  all  are  concentrations 
measured  in  moles  per  liter  (mol/L).  [For  each  additional  proton  a multiprotic  acid  or 
base  can  have,  there  is  similar,  additional  K value.]  The  basis  for  this  definition  is  Le 
Chatelier's  Principle 


[HA]  <=>  [H  + ] + [A~]  (G-2) 

which  states  that  a change  in  the  concentration  of  any  of  the  species  will  cause  a shift  in 

the  equilibrium  to  exactly  preserve  a constant  ratio  (Ka).  The  immediate  consequence  is 

that  each  weak  acid  or  base  has  the  ability  to  drive  the  overall  pH  of  the  system  (since 

they  all  contain  the  term  [H+]),  and  therefore  all  such  equations  are  coupled.  Most 

importantly,  they  are  also  highly  non-linear. 
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APPENDIX  H 

APPLICATION  OF  DESCARTES'  RULE  OF  SIGNS 
To  find  the  nature  of  the  roots  of  the  polynomial  given  in  main  body  of  the  text  as 

(4-24) 


(H*]J+([Na*]  + K,)[H*]!  + 
([Na*]K,-K,[A]-Kw)[Hn-K,Kw=0 
Descartes'  Rule  of  Signs  dictates  examination  of  the  coefficients.  Let 


A = 1 A= 1 

(El-1) 

B = [Na+]  + Ka 

(H-2) 

C = [Na+  ]Ka  - Ka  [A]  - Kw 

(H-3) 

D = -KaKw 

(H-4) 

Due  to  the  fact  that  concentrations  must  always  be  non-negative  and  dissociation 
constants  are  positive,  it  is  concluded  that  A,  B,  and  D always  have  the  same  sign  (+,  +, 
and  -).  Only  C can  vary,  and  negative,  positive,  and  zero  are  all  possibilities.  This 
ensures  that  the  polynomial  will  have  one  of  exactly  three  forms  when  only  the  signs  of 
the  coefficients  are  considered 


X3  +X2  +X-1  = 0 

(H-5) 

X3  + X2-X-1  = 0 

(H-6) 

X3  + X2  -1  = 0 

(H-7) 
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In  essence,  Descartes'  Rule  of  Signs  states  that  there  as  many  positive  real  roots  of 
a polynomial  as  there  are  sign  changes  of  the  coefficients  (Borowski,  1991).  In  addition, 
replacing  X with  (-X)  before  counting  sign  changes  indicates  the  number  of  negative  real 
roots.  The  remaining  roots  must  then  be  complex.  Finally,  zeros  are  ignored. 

Examining  each  of  the  three  cases: 

Case  1 

X3  + X2+X-1  = 0 (H-8) 

There  is  one  sign  change,  indicating  one  positive  real  root;  substituting  (-X) 


-X3+X2-X-1  = 0 (H-9) 

There  are  now  two  sign  changes,  indicating  two  positive  real  roots.  Furthermore,  a cubic 

polynomial  has  exactly  three  roots,  meaning  there  must  be  exactly  one  positive  real  root, 

two  negative  real  roots,  and  zero  complex  roots. 

Case  2 


X3  + X2-X-1  = 0 (H-10) 

Again,  there  is  one  sign  change,  indicating  one  positive  real  root;  substituting  (-X) 


-X3+X2+X-1  = 0 (H-l  1) 

There  are  again  now  two  sign  changes,  indicating  two  positive  real  roots.  Once  more,  it 

is  concluded  that  there  must  be  exactly  one  positive  real  root,  two  negative  real  roots,  and 

zero  complex  roots. 

Case  3 

X3  + X2  -1=0  (H-l  2) 

Again,  there  is  one  sign  change,  indicating  one  positive  real  root;  substituting  (-X) 


-X3+X2 


1 = 0 


(H-13) 
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There  are  again  now  two  sign  changes,  indicating  two  positive  real  roots.  Once  more,  it 
is  concluded  that  there  must  be  exactly  one  positive  real  root,  two  negative  real  roots,  and 
zero  complex  roots. 

All  three  possible  scenarios  result  in  the  same  conclusion — there  is  always  exactly 
one  real  positive  root  and  two  real  negative  roots  for  (4-24). 


APPENDIX  I 

DEVELOPMENT  OF  THE  EXHAUSTIVE  NUMERICAL  SEARCH  USED  TO 
OPTIMIZE  CONTROL  PERFORMANCE 


To  find  the  optimal  controller  gain  Kc  and  integral  time  constant  x,  for  closed- 
loop  performance  with  respect  to  the  ITAE  criterion,  an  adaptive  gridding  system  is  used. 
The  first  grid  range  for  the  controller  gain  can  be  selected  arbitrarily,  and  a decade  of 
potential  solutions  is  used  (e.g.  0.001  to  0.009  by  steps  of  0.001).  For  the  integral  time 
constant,  familiarly  with  the  system  suggests  that  the  decade  of  potential  solutions  from 
10  to  100  by  steps  of  10  should  be  used  for  regulation  responses  and  100  to  1000  by  steps 
of  1 00  for  servo  responses  should  be  used. 

After  the  initial  run  is  completed,  the  controller  gain  Kc  is  determined  to  either 
be  with  the  bounds  or  at  one  of  the  bounds.  If  it  equal  to  one  of  the  bounds,  the  decade  is 
adjusted  by  one  decade  in  the  proper  direction  to  allow  the  potential  solutions  proposed  in 
the  second  run  to  encompass  values  beyond  the  bound  equal  to  the  solution  found  in  the 
first  run  (i.e.,  if  0.001  to  0.009  by  steps  of  0.001  returns  a solution  of  0.009,  then  the  next 
run  should  be  from  0.01  to  0.09  by  steps  of  0.01 ).  This  process  is  repeated  until  the 
optimal  solution  falls  within  the  range  of  bounds  of  the  latest  trial.  For  example,  if  the 
second  iteration  is  from  0.01  to  0.09  by  steps  of  0.01  and  0.03  yields  the  lowest  ITAE,  the 
first  step  in  the  process  is  complete.  In  this  case,  the  optimal  solution  reported  ± the 
previous  step  size  become  the  bounds  for  the  second  step.  For  example,  a solution  of 
0.03  in  steps  of  0.01  will  dictate  bounds  of  0.02  and  0.04  for  the  second  step. 
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The  lone  exception  to  this  scenario  is  when  the  lowest  ITAE  is  found  at  a bound 
of  the  n-th  iteration  and  subsequently  at  the  adjacent  bound  of  the  next  iteration  (e.g. 
when  a range  of  0.001  to  0.009  by  steps  of  0.001  returns  an  optimal  solution  of  0.009  and 
the  next  iteration  returns  an  optimal  solution  of  0.01 ).  If  this  is  the  case,  the  solution  is 
known  to  fall  between  one  appropriate  step  to  the  outside  of  the  two  values,  and  one  step 
outward  from  each  becomes  the  bounds  for  the  second  step  (i.e.,  in  the  scenario 
described,  the  new  bounds  are  0.008  and  0.02). 

With  the  bounds  discovered  for  the  second  step,  the  step  size  is  reduced  to  one 
tenth  of  the  previous  step  size,  and  the  optimization  is  re-run.  For  example,  with  bounds 
of  0.02  and  0.04,  a step  size  of  0.001  is  used.  Similarly,  if  the  bounds  are  determined  by 
the  exception  given  above,  the  step  size  will  be  one  tenth  of  the  finer  of  the  two  bounds 
(i.e.,  in  the  case  of  0.009  (bound  at  0.008)  and  0.01  (bound  at  0.02)  previously  given,  the 
step  size  becomes  0.0001). 

In  this  way,  the  second  step  will  return  a number  for  the  controller  gain  Kc  such 
as  0.23  following  the  first  example.  This  is  said  to  be  accurate  to  two  significant  figures 
as  referred  to  in  the  main  body  of  the  text.  In  the  second  example  given  for  the 
exception,  a plausible  solution  might  be  0.0094. 

With  controller  gain  Kc  identified  to  two  significant  figures,  the  integral  time 
constant  ot  the  controller  xf  is  examined  to  confirm  that  it  is  with  the  chosen  decade  for 
its  possible  solutions.  If  it  is,  the  optimization  is  complete.  If  it  is  not,  the  value  must  be 
100  (recall  that  bounds  of  10  and  1000  are  imposed  on  the  optimization),  and  the  bounds 
are  revised  to  allow  the  optimal  solution  to  include  more  extreme  values  of  x, . During 
this  process,  familiarity  with  the  behavior  of  the  system  dictates  that  the  gain  Kc  should 
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be  expanded  from  the  optimal  solution  by  five  step  sizes  on  either  side  to  ensure  that  fine 
changes  in  the  value  of  controller  gain  Kc  for  the  optimal  solution  are  recorded.  When 
both  parameters  fall  within  the  optimization  bounds,  with  the  gain  Kc  established  to  two 
figures  and  the  time  constant  x,  to  one  figure,  the  optimization  is  complete. 


APPENDIX  J 

MOVING  AVERAGE  AND  EXPONENTIAL  FILTERS 
Moving  Average  Filter  (MAF) 

The  MAF  calculates  the  mean  of  a group  of  previous  data  points.  It  has  one 
parameter,  N,  which  is  the  number  of  points  that  are  included 

yf  (k)  - y(k) + y(k  - 1 ) + y(k  - 2)+-  • -+y(k  - (n  - 1 )) 

y.(k)=TfZy(k-i)  (J-2) 

N i=i 

To  reduce  the  repetitive  calculations,  a more  efficient  version  can  be  derived  by  using  the 
previous  filtered  value 


yf(k-1)  = rrZy(k-i-1)  (J-3) 

N i-0 

Subtracting 


yf(k)-yf(k-l)  = 


N-l 


N-l 


y(k  - 1)  - ^ y(k  - i - 1 ) 


(J-4) 


i=0 


i=0 


yf(k)-yf(k-i)  = 

N-J  N-2  "I  /T 

y(k) + Z y(k  - o - 2 y(k  - 1 - 1 ) - y(k  - (N  - 1 ) - 1 ) 

i=!  i=0 

Shifting  the  index  of  the  second  summation  results  in  the  two  summations  summing  to 
zero,  which  yields 


1 

N 
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yf  (k)  = yf  (k  - 1 ) + ^j-[y(k)  - y(k  - N)]  (J-6) 

Note  that  the  filter  must  be  initialized,  since  the  first  (N-l)  calculations  will  not  have 
enough  data  for  every  term.  In  practice,  this  can  be  done  with  any  reasonable  set  of 
constants,  since  any  error  introduced  from  them  is  quickly  eliminated  by  the  integral 
action  of  the  controller. 

The  implication  is  that  it  requires  a full  N samples  for  a significant  change  in 
input  to  be  completely  realized.  The  response  for  a standard  step  change  is  as  follows 


Figure  J-l.  Illustration  of  how  a step  change  is  observed 
after  passing  through  a moving  average  filter. 

Therefore,  the  delay  manifests  itself  as  a ramp  with  duration  of  N sample  intervals. 
Exponential  Filter  (EF) 

The  EF  calculates  a weighted  average  of  all  previous  data.  Similarly,  it  has  one 
parameter,  a.  It  is  a recursive  relationship,  with  a being  the  weight  given  to  the  new  data 
point,  and  (1-a)  the  weight  assigned  to  the  previously  calculated  value 
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y f (k)  = ay(k)  + (1  -a)yf(k  -1)  (J-7) 

As  new  samples  arrive,  previous  individual  data  points  are  weighted  less  and  less. 

Theoretically,  all  old  data  is  part  of  every  new  calculation,  but  becomes  insignificant  after 
some  time  has  passed 


Figure  J-2.  Illustration  of  how  a step  change  is  observed 
after  passing  through  an  exponential  filter. 

In  contrast  to  the  MAF,  there  is  no  initialization  required,  provided  the  first  sample  is 
simply  passed  through  without  being  averaged. 

As  a result  of  the  algorithm,  the  delay  manifests  itself  in  the  shape  of  an 
exponential  function,  hence  the  name. 
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